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SUMMARY 

This  report  documents  the  area  of  modal  parameter  estimation  in  terms  of  reviewing  efforts  over  the 
past  twenty-five  years,  in  developing  several  new  multiple  reference  methods,  and  in  attempting  to 
provide  a  common  basis  and  understanding  for  all  of  the  modal  parameter  estimation  methods 
developed  to  date.  The  review  of  modal  parameter  estimation  includes  a  substantial  literature  review 
and  the  presentation  of  previous  methods,  such  as  the  Least  Squares  Complex  Exponential,  as  special 
cases  of  general  methods,  such  as  the  Polyreference  Time  Domain  method.  Several  new  modal 
parameter  estimation  methods  are  developed  and  presented  using  consistent  theory  and 
nomenclature.  The  methods  that  are  presented  in  this  manner  include:  Polyreference  Time  Domain, 
Polyreference  Frequency  Domain,  Multiple  Reference  Ibrahim  Time  Domain,  Multiple  Reference 
Orthogonal  Polynomial,  and  Multi  MAC  These  methods,  in  terms  of  general  characteristics,  are  also 
compared  to  other  methods  such  as  the  Least  Squares  Complex  Exponential,  Ibrahim  Time  Domain, 
Eigensystem  Realization  Algorithm,  and  Direct  Parameter  Estimation  methods.  These  methods  are 
all  similar  in  that  the  methods  involve  the  decomposition  of  impulse  response  functions  (time 
domain),  frequency  response  functions  (frequency  domain),  or  forced  response  patterns  (spatial 
domain)  into  characteristic  functions  in  the  appropriate  domain.  These  characteristic  functions  are 
the  single  degree  of  freedom  information  in  the  respective  domain. 
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1.  INTRODUCTION 


Modal  parameter  estimation  is  the  estimation  of  frequency,  damping,  and  modal  coefficients  from  the 
measured  data  which  may  be  in:  (1)  relatively  raw  form  in  terms  of  force  and  response  data  in  the 
time  or  frequency  domain,  or  (2)  in  a  processed  form  such  as  frequency  response  or  impulse  response 
functions.  Most  modal  parameter  estimation  is  based  upon  the  measured  data  being  the  frequency 
response  function;  or  the  equivalent  impulse  response  function,  typically  found  by  inverse  Fourier 
transforming  the  frequency  response  functibn.  Regardless  of  the  form  of  the  measured  data,  the 
modal  parameter  estimation  techniques  have  traditionally  been  divided  into  two  categories:  (1) 
single  degree-of-freedom  (SDOF)  approximations,  and  (2)  multiple  degree-of-freedom  (MDOF) 
approximations.  Since  the  single  degree-of-freedom  equations  are  simply  special  cases  of  the  multiple 
degree-of-freedom  equations,  all  theoretical  discussion  is  made  only  in  terms  of  the  multiple  degree- 
of-freedom  case. 

The  current  effort  in  the  modal  parameter  estimation  area  is  concerned  with  a  unified  theory  that 
explains  any  previously  conceived  modal  parameter  estimation  method  as  a  subset  of  a  general 
theory.  This  unified  theory  concept  would  eliminate  the  confusing  nomenclature  that  currently  exists 
and  simplify  the  understanding  of  the  strengths  and  weaknesses  of  each  method.  The  modal 
parameter  estimation  methods  that  have  been  developed  over  the  past  several  years  involve  multiple 
measurement,  multiple  reference  concepts  that  can  be  viewed  as  an  interaction  between  the  temporal 
domains  (time,  frequency,  etc.)  and  the  spatial  domains  (physical  coordinates,  modal  coordinates, 
etc.)  in  order  to  achieve  the  "best"  estimate  of  the  modal  parameters. 

1.1  Historical  Overview 

While  engineers  have  tried  to  estimate  the  vibration  characteristics  of  structures  since  the  turn  of  the 
centuiy,  the  actual  history  of  experimental  modal  parameter  estimation  is  normally  linked  to  the  work 
by  Kennedy  and  Pancuf1!  in  1947.  Until  this  time,  the  instrumentation  that  had  been  available  was 
not  sufficiently  refined  to  allow  for  detailed  study  of  experimental  modes  of  vibration.  As  the 
instrumentation  and  analysis  equipment  has  improved  over  the  last  forty  years,  major  improvements 
in  modal  parameter  estimation  techniques  have  followed.  Specifically,  the  development  of  accurate 
force  and  response  transducers,  the  development  of  test  equipment  based  upon  digital  computers  and 
the  development  of  the  fast  Fourier  transform  (FFT)  have  been  the  key  advances  that  have  initiated 
bursts  of  development  in  the  area  of  modal  parameter  estimation. 

During  this  time  period,  efforts  in  modal  parameter  estimation  have  involved  two  concepts.  The  first 
concept  involved  techniques  oriented  toward  the  forced  normal  mode  approach  to  modal  parameter 
estimation.  This  approach  to  the  estimation  of  modal  parameters  involves  exciting  the  system  into  a 
single  mode  of  vibration  by  using  a  specific  sinusoidal  forcing  vector.  Since  the  success  of  this 
method  is  determined  by  the  evaluation  of  the  phase  characteristics  with  respect  to  the  characteristics 
occurring  at  resonance,  this  approach  can  be  broadly  classified  as  the  phase  resonance  method.  In 
order  to  refine  the  phase  resonance  method,  particularly  the  force  appropriation  aspect  of  the 
method,  efforts  to  estimate  the  modal  information  on  a  mode  by  mode  basis  using  measured 
impedance,  or  frequency  response,  functions  began.  Much  of  the  early  work  on  this  concept 
centered  on  using  the  phase  information  as  a  means  of  identifying  the  effects  of  separate  modes  of 
vibration  in  the  measurement.  For  this  reason,  this  concept  has  become  known  as  the  phase 
separation  method.  Most  methods  that  are  in  use  today  can  be  classified  as  phase  separation 
methods  since  no  effort  is  made  to  excite  only  one  mode  of  vibration  at  a  time. 

With  respect  to  the  phase  resonance  methods,  the  basic  theory  was  first  documented  by  Lewis  and 
Wrisleyt2)  in  1950.  This  theory  was  refined  and  presented  in  a  more  complete  manner  by  F.  de 
Veubekel3)  in  1956.  Significant  advances  in  the  approach  to  force  appropriation  were  documented  by 
Trail-Nasht4!  and  AsheH5!  in  1958.  Significant  improvements  and  refinements  in  the  phase 


resonance  methods  have  taken  place  in  the  last  thirty  years,  particularly  in  the  automation  of  the 
force  appropriation  and  the  use  of  digital  computers,  but  the  basic  theoretical  concept  has  not 
changed  since  1950. 

With  respect  to  the  phase  separation  methods,  much  effort  has  occurred  over  the  last  forty  years  and 
continues  to  the  present  day.  Kennedy  and  Pancut1]  in  1947  documented  that  the  presence  of  two 
modes  of  vibration  could  be  detected  by  observing  the  rate  of  change  of  the  phase  in  the  area  of 
resonance.  Since  this  method  was  developed  based  upon  a  plot  of  the  real  part  of  the  impedance 
function  versus  the  imaginary  part  of  the  impedance  function,  this  method  is  referred  to  as  a  circle-fit 
method  based  upon  these  characteristics  in  the  Argand  plane.  Broadbentf6!  applied  this  concept  to 
flight  flutter  data  in  1958.  Since  the  data  acquisition  process  was  largely  analog  until  1970,  most  of 
the  work  until  that  time  was  oriented  towards  trying  to  fit  a  single  degree  of  freedom  model  to 
portions  of  the  analog  data.  The  significant  contributors  during  this  time  period  began  with  Stahle  I7! 
in  1958,  and  continued  with  Bishop  and  Giadwell  I8I  and  Pendered  and  Bishop  P-11!  in  1963,  and 
Mahalingham  t12l  in  1967.  Once  data  began  to  be  collected  and  stored  in  a  digital  fashion,  the  phase 
separation  methods  migrated  to  multiple  degrees  of  freedom  approaches.  The  initial  work  involving 
multiple  degree  of  freedom  models  was  documented  by  Klosterman  l13l  in  1971,  Richardson  and 
Potter  I14l  in  1974  and  Van  Loon  1151  in  1974.  While  the  work  during  this  period  involved  the  basic 
polynomial  and  partial  fraction  models  that  are  the  basis  of  modem  experimental  modal  parameter 
estimation  methods,  the  solution  algorithms  were  basically  unstable,  iterative  approaches  to  the 
solution  for  the  unknowns.  Also,  these  methods  used  only  one  measurement  at  a  time  in  the 
estimation  of  the  modal  parameters.  In  1978,  Brown  E16l  documented  work  on  the  Least  Squares 
Complex  Exponential  method  that  was  a  two  stage  approach  to  the  estimation  of  modal  parameters 
using  all  of  the  available  data.  In  the  first  stage,  the  frequency  and  damping  values  are  estimated;  in 
the  second  stage,  the  modal  coefficients  are  estimated.  Ibrahim  I17!,  also  in  1977,  documented  the 
initial  version  of  the  Ibrahim  Time  Domain  Method,  which  formulated  the  solution  for  the  modal 
parameters  into  an  eigenvalue-eigenvector  solution  approach.  These  last  approaches  represent 
conceptual  approaches  that  have  been  extended  today  into  similar  methods  involving  multiple 
references.  "Die  significant  advances  in  the  multiple  reference,  or  polyreference,  methods  used  and 
being  developed  at  the  present  time  were  first  documented  by  Void  l18l  in  1982  with  the 
Polyreference  Time  Domain  method.  Since  that  time,  several  other  polyreference  methods  have 
been  developed.  Detailed  documentation  of  the  multiple  reference  methods  is  contained  in  later 
sections. 

In  summary,  over  the  last  forty  years,  many  experimental  modal  parameter  estimation  methods  have 
been  developed  that  can  be  classified  as  either  phase  resonance  or  phase  separation  methods.  Often, 
it  seems  that  these  methods  are  very  different  and  unique.  In  reality,  the  methods  all  are  derived 
from  the  same  equation  and  are  concerned  with  the  decomposition  of  a  composite  function  into  its 
constituent  parts.  This  decomposition  may  occur  in  the  time  domain  in  terms  of  damped  complex 
exponentials,  in  the  frequency  domain  in  terms  of  single  degree-of-freedom  functions,  or  in  the 
modal  domain  in  terms  of  modal  vectors.  This  decomposition  may  occur  during  the  test,  as  in  the 
phase  resonance  methods,  or  occur  during  analysis,  as  in  the  phase  separation  methods.  The  various 
modal  parameter  estimation  methods  are  enumerated  in  the  following  list: 

•  Forced  Normal  Mode  Method  t2'5-19-20) 

•  Quadrature  Amplitude  t7*8,11! 

.  Kennedy-Pancu  Circle  Fit  H.D, 21-23] 

.  Single  Degree-of-Freedom  Polynomial  l14-21-22-24-25) 

•  Nonlinear  Frequency  Domain  I13’14-21’22! 

•  Complex  Exponential  ^26,27^ 


•  Least-Squares  Complex  Exponential  (LSCE) 

•  Ibrahim  Time  Domain  (11D)  l28"32^ 

•  Eigensystem  Realization  Algorithm  (33'3S1 

•  Orthogonal  Polynomial  ^36,37^ 

•  Global  Orthogonal  Polynomial 

•  Polyreference  Time  Domain  t18*39*40! 

•  Polyreference  Frequency  Domain  l40-43! 

•  Direct  Parameter  Identification:  Time  Domain  ^ 

•  Autoregressive  Moving  Average  (ARMA)  t44"48! 

•  Direct  Parameter  Identification:  Frequency  Domain  ^40,49^ 

12  Multiple-Reference  Terminology 
12.1  Mathematical  Models 

The  most  general  model  that  can  be  used  is  one  in  which  the  elements  of  the  mass,  damping,  and 
stiffness  matrices  are  estimated,  based  upon  measured  forces  and  responses.  Thus,  the  model  that  is 
used  is  based  upon  a  matrix  differential  equation  transformed  into  the  domain  of  interest. 

Time  domain: 

{*}♦[*]  {*}={/}  (1) 

Frequency  domain: 

-JIM)  W  *  MC]  {X)  ♦  [K]  { X }  =  {F}  (2) 

Laplace  domain: 

sa[M]  {X}  +  s[C]  {X}  ♦  [X]  W  =  {F}  (3) 

If  Eq.  (1),  (2),  or  (3)  is  used  as  the  model  for  parameter  estimation,  the  elements  of  the  unknown 
matrices  must  first  be  estimated  from  the  known  force  and  response  data  measured  in  the  time  or 
frequency  domain.  Once  the  matrices  have  been  estimated,  the  modal  parameters  can  be  found  by 
the  solution  of  the  classic  eigenvalue-eigenvector  problem  ^42,49l  Due  to  truncation  of  the  data  in 
terms  of  frequency  content,  limited  numbers  of  degrees -of-freed om,  and  measurement  errors,  the 
matrices  found  by  Eq.  (1),  (2),  or  (3)  are,  in  general,  not  directly  comparable  to  matrices  determined 
from  a  finite  element  approach.  Instead,  the  matrices  that  are  estimated  simply  yield  valid  input- 
output  relationships  and  valid  modal  parameters.  This  is  because  there  is  an  infinite  number  of  sets 
of  mass,  damping,  and  stiffriess  matrices  that  yield  the  same  modal  parameters  over  a  reduced 
frequency  range  limited  to  the  dynamic  range  of  the  measurements.  For  this  reason,  Eqs.  (1),  (2),  and 
(3)  are  often  pre-multipied  by  the  inverse  of  the  mass  matrix  so  that  the  elements  of  the  two  matrices 
[D]  and  [E]  are  estimated: 


Time  domain: 


[/]{*W0]{iM£H*}«{/'>  (4) 

Frequency  domain: 

+?[I]{X}*MD]{X}  +  [E]{X}~{F'}  (5) 

Laplace  domain : 

s*[I]{X}+s[D]{X}  +  [E]{X}-{F'}  (6) 

Existing  modal  parameter  estimation  methods  used  in  commercial  modal  analysis  systems  most  often 
employ  a  model  based  upon  measured  impulse  response  (time  domain)  or  frequency  response 
(frequency  domain)  functions.  While  the  exact  model  used  as  the  basis  for  modal  parameter 
estimation  varies,  almost  ail  models  used  in  conjunction  with  frequency  response  function  data  can  be 
described  by  a  general  model  in  the  time  domain,  frequency  domain,  or  Laplace  domain.  The  general 
model  in  the  time  domain  is  a  damped  complex  exponential  model  (often  the  impulse  response 
function)  while  the  general  model  in  the  frequency  domain  is  the  frequency  response  function.  The 
general  model  in  the  Laplace  domain  is  the  transfer  function.  For  general  viscous  damping,  the 
mathematical  models  for  each  domain  for  a  multiple  degree  degree-of-freedom  mechanical  system 
can  be  stated  as: 

Time  Domain: 

(7) 

r=d 


where: 

s  =  Laplace  variable 

s  =  a  ♦  ju 

o  —  angular  damping  variable  (rad/sec) 
w  —  angular  frequency  variable  (rad/sec) 
p  =  measured  degree-of-freedom  (response) 

q  —  measured  degree-of-freedom  (input) 

r  —  modal  vector  number 

N  —  number  of  modal  frequencies 

—  residue 
Apr-  “  QrtMv 

Q,  —  complex  modal  scaling  coefficient  for  mode  r 
i>f,  =  modal  coefficient  for  measured 


degree-of-freedom  p  and  mode  r 
Lr  —  modal  participation  factor  for  reference 
degree-of-freedom  q  and  mode  r 
A,  =  system  pole 
A,  "  o,+M 


The  models  described  in  Eqs.  (7)  through  (9)  have  many  other  equivalent  forms  based  upon 
expansion  of  the  terms  under  the  summation.  Also,  the  models  take  on  slightly  different  forms  under 
assumptions  concerning  specific  physical  damping  mechanisms  (hysteretk,  etc)  [  Other 

forms  of  these  models  are  also  used  where  certain  assumptions  or  mathematical  relationships  are 
utilized.  For  example,  an  equivalent  model  can  be  found  when  the  common  denominator  of  Eq.  (8) 
is  formed  yielding  a  polynomial  numerator  and  polynomial  denominator  of  maximum  order  "2N" 
[13,14,22]  jjjg  denominator  polynomial  is  then  a  function  of  the  system  poles.  Often,  an  assumption  is 
made  concerning  the  modal  vectors  being  normal  (real)  rather  than  complex.  This  reduces  the 
number  of  unknowns  that  must  be  estimated  by  "N". 

122  Sampled  Data 

The  mathematical  models  described  in  the  previous  section  are  all  developed  based  upon  the  concept 
that  the  data  is  continuous.  In  reality  the  data  that  is  available  must  be  thought  of  as  sampled  data  in 
each  domain.  This  restriction  requires  special  consideration  when  applying  the  models  developed  in 
Eqs.  (1)  through  (9).  Differential  equations  must  now  be  thought  of  as  finite  difference  equations; 
continuous  integral  transforms  are  replaced  by  discrete  transforms  such  as  the  Fast  Fourier 
Transform  (FFT)  and  the  Z  Transform.  The  concepts  affecting  the  numerical  processing  of  sampled 
data  with  respect  to  the  continuous  models  represented  in  Eqs.  (1)  through  (9)  are  exactly  the  same 
as  the  concepts  that  are  the  basis  of  the  area  of  digital  signal  analysis  with  respect  to  the 
measurement  of  the  data.  The  limitations  of  the  frequency  information  creates  special  processing 
problems  that  are  related  to  Shannon’s  Sampling  Theorem;  the  limitations  of  the  dynamic  range  of 
the  measured  data  and  of  the  computer  precision  yield  special  numerical  problems  with  respect  to  the 
solution  algorithm. 

In  general,  the  numerical  considerations  often  determine  which  mathematical  model  will  be  most 
effective  in  the  estimation  of  modal  parameters.  Time  domain  models  tend  to  provide  the  best 
results  when  a  large  frequency  range  or  large  numbers  of  modes  exist  in  the  data.  Frequency  domain 
models  tend  to  provide  the  best  results  when  the  frequency  range  of  interest  is  limited  and  when  the 
number  of  modes  is  small.  While  these  are  general  considerations,  the  actual  numerical 
implementation  determines  the  ability  of  the  algorithm  to  estimate  modal  parameters  accurately  and 
efficiently. 

1.2.3  Consistent  Data 

Modal  parameter  estimation  methods  all  assume  that  the  system  that  is  being  investigated  is  linear 
and  time  invariant.  While  this  is  often  nearly  true,  these  assumptions  are  never  exactly  true. 
Consistent  data  refers  to  the  situation  where  the  data  is  acquired  so  as  to  best  satisfy  these  two 
assumptions.  Problems  associated  with  linearity  can  be  minimized  by  maintaining  a  prescribed  force 
level  and/or  using  excitation  methods  that  give  the  best  linear  approximation  to  the  nonlinear 
characteristic  (random  excitation).  Problems  associated  with  the  time  invariance  constraint  can  be 
minimized  by  acquiring  all  of  the  data  simultaneously  using  multiple  excitations  ‘  \  This  reduces 

mass  loading  and  boundary  condition  variations  that  can  be  caused  by  moving  a  transducer  around 
the  structure  or  by  changing  the  location  of  the  excitation. 


1.2.4  Residuals 


Continuous  systems  have  an  infinite  number  of  degrees-of-ffeedom  but.  in  general,  only  a  finite 
number  of  modes  can  be  used  to  describe  the  dynamic  behavior  of  a  system.  The  theoretical  number 
of  degrees-of-ffeedom  can  be  reduced  by  using  a  finite  frequency  range  (/«,/»).  Therefore,  for 
example,  the  frequency  response  function  can  be  broken  up  into  three  partial  sums,  each  covering  the 
modal  contribution  corresponding  to  modes  located  in  the  frequency  ranges  (0,/#),  (/.,/»),  and 
(fk, oo)  as  shown  in  Figure  1. 


FREQUENCY,  HZ 


Figure  1.  Frequency  Range  of  Interest 


In  the  frequency  range  of  interest,  the  modal  parameters  can  be  estimated  to  be  consistent  with  Eq. 
(8).  In  the  lower  and  higher  frequency  ranges,  residual  terms  can  be  included  to  handle  modes  in 
these  ranges.  In  this  case,  the  general  frequency  response  function  model  can  be  stated: 


A 


(10) 


♦*/>« 


where: 

—  residual  effect  of  lower  frequency  modes 
~  residual  effect  of  higher  frequency  modes  (constant  with  w) 


In  many  cases  the  lower  residual  is  called  the  inertia  restraint,  or  residual  inertia,  and  the  upper 
residual  is  called  the  residual  flexibility  (13l  In  this  common  formulation  of  residuals,  both  terms  are 
real-valued  quantities.  The  lower  residual  is  a  term  reflecting  the  inertia  or  mass  of  the  lower  modes 
and  is  an  inverse  function  of  the  frequency  squared.  The  upper  residual  is  a  term  reflecting  the  the 
flexibility  of  the  upper  modes  and  is  constant  with  frequency.  Therefore,  the  form  of  the  residual  is 
based  upon  a  physical  concept  of  how  the  system  poles  below  and  above  the  frequency  range  of 
interest  affects  the  data  in  the  range  of  interest  As  the  system  poles  below  and  above  the  range  of 
interest  are  located  in  the  proximity  of  the  boundaries  of  the  frequency  range  of  interest,  these  effects 
are  not  the  simple  real-valued  quantities  noted  in  Eq.  (10).  In  these  cases,  residual  modes  may  be 
included  in  the  model  to  partially  account  for  these  effects.  When  this  is  done,  the  modal  parameters 
that  are  associated  with  these  residual  poles  have  no  physical  significance,  but  may  be  required  in 
order  to  compensate  for  strong  dynamic  influences  from  outside  the  frequency  range  of  interest. 
Using  the  same  argument,  the  lower  and  upper  residuals  can  take  on  any  mathematical  form  that  is 
convenient  as  long  as  the  lack  of  physical  significance  is  understood.  Power  functions  of  frequency 
(zero,  first,  and  second  order)  are  commonly  used  within  such  a  limitation.  In  general,  the  use  of 
residuals  is  confined  to  frequency  response  function  models.  This  is  primarily  due  to  the  difficulty  of 
formulating  a  reasonable  mathematical  model  and  solution  procedure  in  the  time  domain  for  the 
general  case  that  includes  residuals. 


i 


J 


1.2.5  Global  Modal  Parameters 


Theoretically,  modal  parameters  are  considered  to  be  unique  based  upon  the  assumption  that  the 
system  is  linear  and  time  invariant.  Therefore,  the  modal  frequencies  can  be  determined  from  any 
measurement  and  the  modal  vectors  can  be  determined  from  any  reference  condition.  If  multiple 
measurements  or  reference  conditions  are  utilized,  the  possibility  of  several,  slightly  different,  answers 
for  each  modal  parameter  exists.  The  concept  of  global  modal  parameters,  as  it  applies  to  modal 
parameter  estimation,  means  that  there  is  only  one  answer  for  each  modal  parameter  and  that  the 
modal  parameter  estimation  solution  procedure  enforces  this  constraint.  Every  frequency  response 
or  impulse  response  function  measurement  theoretically  contains  the  information  that  is  represented 
by  the  characteristic  equation  (modal  frequencies  and  damping).  If  individual  measurements  are 
treated  in  the  solution  procedure  independent  of  one  another,  there  is  no  guarantee  that  a  single  set 
of  modal  frequencies  and  damping  are  generated.  In  a  like  manner,  if  more  than  one  reference  is 
measured  in  the  data  set,  redundant  estimates  of  the  modal  vectors  can  be  estimated  unless  the 
solution  procedure  utilizes  all  references  in  the  estimation  process  simultaneously.  Most  modal 
parameter  estimation  algorithms  estimate  the  modal  frequencies  and  damping  in  a  global  sense  but 
few  estimate  the  modal  vectors  in  a  global  sense. 


1.2.6  Modal  Participation  Factors 


A  modal  participation  factor  is  a  complex-valued  scale  factor  that  is  the  ratio  of  the  modal  coefficient 
at  one  reference  degree -of-freedom  to  the  modal  coefficient  at  another  reference  degree-of-freedom. 
A  more  general  view  of  the  modal  participation  factor  is  that  it  represents  the  relationship  between 
the  residue  and  the  eigenvector  coefficient  as  in  the  following  equations: 


-7- 
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(13) 

where: 

p  —  measured  degree-of-freedom  (response) 

q  —  measured  degree-of-freedom  (reference) 

r  —  modal  vector  number 

—  residue 

Qr  —  complex  modal  scaling  coefficient  for  mode  r 

Vy  —  modal  coefficient  for  measured 

degree-of-freedom  p  and  mode  r 

Vy  =  modal  coefficient  for  reference 

degree-of-freedom  q  and  mode  r 

Lr  —  modal  participation  factor  for  reference 
degree-of-freedom  q  and  mode  r 

From  a  mathematical  standpoint,  the  modal  participation  matrix  is  equal  to  the  left  eigenvectors  of 
next  equation : 

M-[*]  r*jH 

where: 

|//j  —  transfer  function  matrix 

|  ¥  j  =  complex  modal  vector  matrix 
[  AJ  =  diagonal  matrix  with  poles 


Note  that  for  Eq.  (12)  the  modal  participation  factor  represents  the  product  of  a  modal  scaling 
coefficient  and  another  term  from  the  right  eigenvector  for  reference  degree-of-freedom  q.  This  will 
always  be  true  for  reciprocal  systems  since  the  left  and  right  eigenvectors  for  a  given  mode  are  equal. 
For  non-reciprocal  systems,  the  modal  participation  factor  is  the  appropriate  term  from  the  left 
eigenvector.  Note  also  that  the  modal  participation  factor,  since  it  is  related  to  the  eigenvector,  has 
no  absolute  value  but  is  relative  to  the  magnitudes  of  the  other  elements  in  the  eigenvector. 

Modal  participation  factors  reflect  the  interaction  of  the  spatial  domain  with  the  temporal  domain 
(time  or  frequency).  Modal  participation  factors  can  be  computed  anytime  that  multiple  reference 
data  are  measured  and  such  factors  are  used  in  multiple  reference  modal  parameter  estimation 
algorithms.  Modal  participation  factors  relate  how  well  each  modal  vector  is  excited  from  each  of  the 
reference  locations.  This  information  is  often  used  in  a  weighted  least  squares  error  solution 
procedure  to  estimate  the  modal  vectors  in  the  presence  of  multiple  references.  Theoretically,  these 
modal  participation  factors  should  be  in  proportion  to  the  modal  coefficients  of  the  reference  degrees 
of  freedom  for  each  modal  vector.  Modal  participation  factors  in  a  solution  procedure  enforce  the 
constraint  concerned  with  Maxwell’s  reciprocity  between  the  reference  degrees  of  freedom.  Most 
multiple  reference,  modal  parameter  estimation  methods  estimate  modal  participation  factors  as  part 
of  the  first  stage  estimation  of  global  modal  frequencies  and  damping. 


1.2.7  Order  of  the  Model 


The  order  of  the  model  equals  the  number  of  unknowns  that  must  be  estimated  in  the  model.  In  the 
modal  parameter  estimation  case,  this  refers  to  the  frequency,  damping,  and  complex  modal 
coefficient  for  each  mode  of  vibration  at  eveiy  measurement  degree-of-freedom  plus  any  residual 
terms  that  must  be  estimated.  Therefore,  the  order  of  the  model  is  directly  dependent  on  the  number 
of  modal  frequencies,  "N”,  that  are  to  be  estimated.  For  example,  for  a  system  with  "N”  modes  of 
vibration,  assuming  that  no  residuals  were  required,  M4NM  unknowns  must  be  estimated.  For  cases 
involving  measured  data,  the  order  of  the  model  is  extremely  important.  Estimates  of  modal 
parameters  are  affected  by  the  order  of  the  model.  A  problem  arises  from  the  inability  to  be  certain 
that  the  correct  order  of  the  model  has  been  chosen  during  the  initial  estimation  phase.  If  the 
number  of  modes  of  vibration  is  more  or  less  than  MN",  modes  of  vibration  will  be  found  that  do  not 
exist  physically  or  modes  of  vibration  will  be  missed  that  actually  do  occur.  In  addition,  the  values  of 
frequency,  damping,  and  complex  modal  coefficient  for  the  actual  modes  of  vibration  will  be  affected. 

The  number  of  modes  of  vibration  is  normally  chosen  between  one  and  an  upper  limit,  dependent  on 
the  memory  limitations  of  the  computational  hardware.  The  true  number  of  system  poles  is  a 
function  of  the  frequency  range  of  the  measurements  used  to  estimate  the  modal  parameters.  By 
observing  the  number  of  peaks  in  the  frequency  response  function,  the  minimum  number  of  system 
poles  can  be  estimated.  This  estimate  is  normally  low,  based  upon  poles  occuring  at  nearly  the  same 
frequency  (psuedo-repeated  roots),  limits  on  dynamic  range,  and  poorly  excited  modes.  For  these 
reasons,  the  estimate  of  the  correct  order  of  the  model  is  often  in  error.  When  the  order  of  the  model 
is  other  than  optimum,  the  estimate  of  the  modal  parameters  will  be  in  error. 

Many  of  the  parameter  estimation  techniques  that  are  used  assume  that  only  one  mode  exists  in  a 
limited  range  of  interest  and  all  of  the  other  modes  appear  as  residual  terms.  For  this  case  Eq.  (10) 
can  be  rewritten  as: 
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U.8  Solution  Procedure 

Equations  (7)  through  (9)  represent  a  nonlinear  model  in  terms  of  the  unknown  modal  parameters. 
This  can  be  noted  from  the  unknowns  in  the  numerator  and  denominator  of  Eq.  (8)  and  the 
unknowns  as  the  argument  of  the  transcendental  functions  of  Eq.  (7).  The  nonlinear  aspect  of  the 
model  must  be  treated  in  one  of  two  ways:  (1)  by  the  use  of  an  iterative  solution  procedure  to  solve 
the  nonlinear  estimation  problem,  allowing  all  modal  parameters  to  vary  according  to  a  constraint 
relationship  until  an  error  criterion  reaches  an  acceptably  low  value,  or  (2)  by  separating  the 
nonlinear  estimation  problem  into  two  linear  estimation  problems.  For  the  case  of  structural 
dynamics,  the  common  technique  is  to  estimate  "2N"  frequencies  and  damping  values  in  a  first  stage 
and  then  to  estimate  the  H4N"  modal  coefficients  plus  any  residuals  in  a  second  stage. 

In  the  iterative  technique  in  the  solution  of  the  nonlinear  estimation  approach,  a  set  of  starting  values 
must  be  chosen  to  initiate  the  sequence.  The  number  and  value  of  these  starting  values  affect  the  final 
result.  Poor  initial  estimates  can  lead  to  problems  of  convergence,  as  a  result  of  which,  close  operator 
supervision  usually  is  required  for  a  successful  use  of  this  technique. 

An  alternative  method  is  to  reformulate  the  nonlinear  problem  into  a  number  of  linear  stages  so  that 
each  stage  is  stable.  The  actual  data  that  are  used  in  the  estimate  of  the  modal  parameters  also  affect 
the  results.  Based  on  the  choice  of  the  order  of  the  model,  "N”,  there  are  "4N"  modal  parameters  to 
be  estimated.  If  residuals,  in  one  form  or  another,  are  also  included,  the  number  of  modal 
parameters  to  be  estimated  will  be  slightly  higher.  The  common  method  of  solving  for  these 


unknown  modal  parameters  is  to  find  an  equation  involving  known  information  for  every  unknown 
that  will  be  found.  In  this  case,  the  frequency  response  function  or  impulse  response  function  that 
has  been  measured  provides  the  known  information  and  Eq.  (7)  or  Eq.  (8)  can  be  repeated  for 
different  frequencies  or  time  values  in  order  to  obtain  the  sufficient  number  of  equations.  These 
equations,  for  the  linear  case,  can  then  be  solved  simultaneously  for  the  unknown  modal  parameters. 
As  an  illustration  of  this  relationship,  consider  a  common  modal  parameter  situation  in  which  the 
number  of  modes  in  the  frequency  range  of  interest  is  between  1  and  30.  Assuming  the  highest 
ordered  model  means  that  slightly  more  than  120  modal  parameters  must  be  estimated.  From  a  single 
frequency  response  measurement,  1024  known  values  of  the  function  will  be  available  (S12  complex 
values  at  succesive  values  of  frequency).  Many  more  equations,  based  on  the  known  values  of 
frequency  response,  can  be  formed  than  are  needed  to  find  the  unknown  modal  parameters.  An 
obvious  solution  is  to  choose  enough  equations  to  solve  for  the  modal  parameters.  The  problem 
arises  in  determining  what  part  of  the  known  information  is  to  be  involved  in  the  solution.  As 
different  portions  of  the  known  data  (data  near  a  resonance  compared  to  an  anti-resonance,  for 
example)  are  used  in  the  solution,  the  estimates  of  modal  parameters  vary.  As  the  quality  of  the  data 
becomes  marginal,  this  variance  can  be  quite  large.  When  the  modal  parameters  that  are  estimated 
appear  to  be  non-physical,  this  is  often  the  reason. 

To  solve  this  problem,  all  or  a  large  portion  of  the  data  can  be  used  if  a  pseudo-inverse  type  of 
solution  procedure  is  used.  One  procedure  that  is  used  is  to  formulate  the  problem  so  as  to  minimize 
the  squared  error  between  the  data  and  the  estimated  model.  This  least-squares  error  method  to  the 
solution  is  the  most  commonly  used  technique  in  the  area  of  modal  parameter  estimation.  If  there 
are  many  more  known  pieces  of  information  than  unknowns  that  must  be  estimated,  many  more 
equations  can  be  formed  than  are  needed  to  solve  for  the  unknowns.  The  least-squares  error  method 
to  the  solution  allows  for  all  of  these  redundant  data  to  be  used  to  estimate  the  modal  parameters  in 
a  computationally  efficient  manner.  The  least-squares  error  method  usually  can  be  derived  directly 
from  the  linear  equations  using  a  normal  equations  approach.  In  genera],  this  procedure  does  not 
increase  the  memory  or  computational  requirements  of  the  computational  hardware  significantly. 
Any  solution  procedure  that  can  be  used  is  only  estimating  a  "best”  solution  based  upon  the  choice  of 
the  model,  the  order  of  the  model,  and  the  known,  measured  data  used  in  the  model. 

1.2.9  Characteristic  Polynomial 

The  impetus  of  this  section  is  to  show  that  for  discrete  data,  a  difference  characteristic  equation  can 
be  formulated  in  order  to  solve  for  the  poles  of  the  system.  Further,  it  will  be  shown  that  the 
difference  equation  can  be  formulated  directly  from  the  impulse  response  function  data.  By  solving 
for  the  polynomial  coefficients  and  the  roots  of  the  polynomial  equation,  the  modal  parameters, 
frequency  and  damping,  are  determined.  The  characteristic  polynomial  will  be  formulated  for  the 
continuous  case,  as  a  differential  function,  and  then  extended  to  the  discrete  case,  as  a  difference 
function. 

1.2.9.1  Differential  Theory 

The  homogeneous  differential  equation  for  a  single  degree  of  freedom  system  is: 

mx(t)  +  cx(t)  +kx(t)  *  0  .  (15) 

In  order  to  solve  the  differential  equation  assume  a  solution  of  the  of  the  form  x  (t)  =  X  e'  *,  where  X 
is  a  scalar  value.  Substituting  the  appropriate  derivatives  of  the  assumed  solution  into  Eq.(15): 

(msa+cs+k)Xe',  =  0 

Thus,  the  differential  equation  is  transformed  into  an  algebraic  polynomial  equation,  called  the 


CTBgggamm gg aftnurnmn v '>"-•  vw v  *  ? r. «r. 


characteristic  equation. 

mj’+n+t*  0  (16) 

The  complex  valued  roots  of  the  characteristic  equation  will  yield  the  characteristic  solutions,  Aj  and 
Aa.  The  real  part  is  the  damping  and  the  imaginary  part  is  the  eigenfrequency,  or  damped  natural 
frequency.  Thus,  the  solution  to  the  governing  differential  equation  is: 

*(0  =  E*  e*1*  . 


The  scalar  magnitudes,  Xx  and  X3,  are  determined  from  the  system  initial  conditions.  Note  that  any 
exponential  function  will  satisfy  the  differential  equation.  One  such  function  of  particular  interest,  is 
the  impulse  response  function. 

*<0-  • 

r= 1 

A  system  with  TV  degrees-of-freedom  can  be  described  by  a  set  of  TV,  coupled,  second  order 
differential  equations.  The  characteristic  equation  for  this  system  is  represented  by  the  following 
polynomial: 


sw  ♦  aw. i  s"'1  ♦  a-MJx  s®-5  s  +  a0  *  0  (17) 

Solving  this  polynomial  equation  will  yield  IN  complex  valued  roots,  or,  characteristic  solutions  Ar. 
Then  the  solutions  to  the  differential  equations  will  be  complex  exponentials  of  the  form: 

2N 

r=l 

where: 

•  p  -  response  location  degree-of-freedom 

•  q  -  reference  location  degree-of-freedom. 

Thus,  impulse  response  functions, 

MO  =  EM  *A,‘ 

r=l 


1 


will  also  satisfy  the  differential  equations.  Consider  a  few  impulse  response  functions  for  various 
reference  and  response  points. 

*u(0-  EM  <?Ar‘ 

r* I 


h  12(0«  EM 

r*l 

*1 3(0*  EM  ^ 

r=l 


mo- 


2N 


E 


t 


X,t 
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The  common  characteristic  in  each  of  the  above  equations  is  that  every  impulse  response  function  is  a 
linear  superposition  of  identical  damped  complex  exponentials,  ex,‘  for  /■  »  1  ->  2N.  That  is,  the  roots 
of  the  characteristic  polynomial  are  common  to  all  reference  and  response  locations.  Thus,  the 
characteristic  solutions  are  global  system  parameters,  since  they  are  independent  of  reference  or, 
response  location.  The  important  result  is  that  since  each  ex,i  is  a  characteristic  solution  to  the 
homogeneous  linear  differential  equation, 

* aw.7D^[A^eKt]  ♦. ..  d[a„  eM)  +  ao  =  0 

where, 

•  D*[f  (01  =  (differential  operator) 

dt “ 

that  a  linear  superposition  of  characteristic  solutions  will  also  be  a  solution.  That  is,  hm(t)  will  also 
satisfy  the  differential  equation.  Actually,  a  set  of  N  second  order  linear  differential  equations  must 
be  satisfied,  but,  a  differential  equation  of  order  IN  can  be  found  that  will  have  the  same  roots  as  the 
set  of  N  second  order  equations. 

Substituting  a  few  impulse  response  functions  for  various  reference  and  response  points,  a  number  of 
differential  equations  are  obtained. 


a2S  +  a2N-l  D2*'1  ^Au(/)j  + 

♦«,  +«o  =  0 

a?nDw^h u(f)j  *aM r.i  D2*'1  ^/»ia(0]  ♦  fla*r-a  ♦  ■ 

.  •  +ai  D^/iia(f)]  +a0  =  0 

a  aw  +<*aw.i  ♦flaw-a  +• 

■  •  ♦<*!  £>(Ais(0]  ♦  a0  =  0 

awDw[h„(t)}  ♦  flaw-iO3*’1  (*„(')]  ♦%.aOM-,(Ap,(0) 

•  ♦«!  *a0  =  0 

Note  that  the  coefficients  o0toaudo  not  vary  with  reference  or  response  location  and  thus,  can  be 
estimated  from  a  combination  of  various  number  of  reference  and  response  points. 

1.2.9 .2  Difference  Theory 


From  an  experimental  standpoint,  the  data  are  sampled,  which  means  instead  of  continuous 
knowledge  of  the  system,  the  values  obtained  are  for  distinct  discrete  time  points.  The  impulse 
response  functions  are  generally  obtained  by  inverse  Fourier  transforming  the  frequency  response 
functions.  Thus,  from  the  discrete  impulse  response  functions  the  pole  information,  frequency  and 
damping,  is  determined.  The  model  for  the  discrete  impulse  response  function  is: 

*„('*)  =  E'V  cMi  =  E Apr*  (18) 

rsl  rsl 


where, 

•  /*  »  k  At 

•  A  /  is  the  sample  interval 

•  blocksize 
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It  should  be  noted,  for  discrete  data,  that  the  sample  interval,  A  t,  limits  the  frequency  for  which  valid 
information  can  be  determined,  whereas,  in  the  analysis  of  continuous  data,  there  are  no  frequency 
constraints.  In  other  words,  theoretically,  characteristics  can  be  determined  to  infinite  frequency  for 
continuous  functions,  but,  the  process  of  digitally  sampling  continuous  data  causes  a  maximum 
frequency  for  which  characteristics  can  be  determined.  The  frequencies  above  this  maximum  will 
alias  back  into  the  sampled  bandwidth,  and  thus  bias  the  results.  For  this  reason,  low-pass  filters  are 
used  to  exclude  information  above  the  maximum  frequency. 


Recall  the  characteristic  equation  for  an  n  degree-of-freedom  system: 

flaw-i  +  flaw  s™  m  0 


(19) 


will  have  2 N  characteristic  solutions,  A„  for  r  =  1  — ►  2N.  The  characteristic  polynomial  is  not  unique 
in  that,  many  polynomials  can  be  constructed  that  will  yield  the  same  characteristic  solutions,  even 
though  the  coefficients  will  be  different  For  this  reason,  another  polynomial  can  be  formulated  that 
will  have  characteristic  solutions  that  are  related  to  the  characteristic  solutions  of  Eq.(19).  The 
polynomial  has  the  form: 

a'0  z°  ♦  a\  z1  ♦  a\  z3  ■*■...  +  a^y.i  z3**1  ♦  dw  z w  =  0  .  (20) 

The  relationship  between  z  and  s  is  z  *  e*A‘.  Analogous  to  Eq.(19),  there  are  also  2N  characteristic 
solutions  of  Eq.(20),  z,  for  r  =  1  -*  2N.  The  roots  of  the  two  equations  are  related  by  z,  =  eKCkt 
where,  zr  are  precisely  the  values  of  z  for  which  the  characteristic  equation,  Eq.(20),  is  zero.  Note  that 
z,  is  simply  the  sampled  form  of  the  continuous  exponential  solution  in  the  differential  case.  Thus,  by 
knowing  the  system  characteristics,  z,,  the  desired  parameters,  Ar,  can  be  determined.  If  the 
coefficients  are  known,  Eq.(20)  could  be  solved,  but,  from  an  experimental  aspect,  both  the 
coefficients  and  the  system  characteristics  are  unknown.  Thus,  in  order  to  determine  the  system 
characteristics,  the  a  coefficients  must  be  determined  first.  This  is  accomplished  by  substituting  a 
characteristic  solution  of  the  system,  z,,  into  Eq.(20). 


flo  z°  +  a\  zi  *  a\  A„  z3  ♦ . . .  ♦  dw.x  Anr  z™-1 
Substituting  z,  =  ex,At  into  the  above  equation, 

'  j  /  At  \0  ,  A  /  \1  .  '  a  /  ArAt  . 


*  a1S  /*p*r  -  0 

...♦%V(^')” 


=  0 


or, 


Aptre°  +  a\  A. 


.ex’&t*d3Af 


.eA'3AV 


*  aw^pr  e 


X,2SAt  _ 


(21) 


The  important  result  is  that  since  each  eA'A‘  is  a  characteristic  solution  to  the  homogeneous  linear 
difference  equation,  that  a  linear  superposition  of  characteristic  solutions  will  also  be  a  solution  of 
Eq.(21),  which  means  that,  in  general,  Eq.(18)  can  be  substituted  into  Eq.(21).  Once  again,  a  set  of  N 
second  order  linear  difference  equations  must  be  satisfied,  but,  a  difference  equation  of  order  IN  can 
be  found  that  will  have  the  same  roots  as  the  set  of  N  second  order  equations. 


Consider  a  number  of  equations  for  various  reference  and  response  locations: 


A 


A 


*1 


a 


i 


1 


si 


ao  ^ufro)  +  ai  hn(ti)  +  a3  hu(t a)  * . . .  +  flay. i  /tn(f2w- i)  ♦  flaw  ^n(^av)  =  0 


ao  ^iafro)  *  ai  ^ia(fi)  *  ai  h  i3(,t3)  ♦. . .  ♦  flaw-i  ^ia(taw-i)  *  flaw  ^ia(faw)  =  0 
ao  ^is(to)  +  ai  +  fla  hl3(t2)  •  •  ♦  flaw-i  ^is(faw-i)  +  aw  hl3(t2 y)  =  0 

ao  hpfito)  *  ai  hpf(ti)  ♦  a3  hpfO 2)  + . . .  ♦  flay-i  fy*  (faw-i)  *  flaw  (taw)  =  0 


i 
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Note  that  the  coefficients  a'0  to  a\n  do  not  vary  with  reference  or  response  location  and  thus,  can  be 
estimated  from  a  combination  of  various  number  of  reference  and  response  points.  Once  the  a 
coefficients  are  estimated  from  a  set  of  equations  similar  to  the  ones  above,  the  poles,  zr,  and  hence 
Ar,  can  be  estimated  from  the  2N  solutions  of  the  characteristic  equation, 

a'0  z°  +  a\  z1  +  a2  z2  + . . .  +  aw.x  z +  a'w  zw  =  0 

where: 


In  summary,  a  series  of  2N  linear  difference  equations  with  constant  coefficients  are  formed  from  the 
sampled  impulse  response  function  data  in  order  to  solve  for  the  common  constant  coefficients. 
These  coefficients  are  then  used  in  the  characteristic  equation  to  solve  for  the  system  characteristics, 
zT,  which  contain  the  desired  parameters,  Ar. 

Note  that  the  characteristic  polynomial  for  the  continuous,  or  discrete  case,  is  of  order  2 N,  that  is, 
twice  the  number  of  modes.  This  results  in  a  time  domain  differential,  or  difference  equation  of 
ovder  2N.  For  this  reason,  from  a  numerical  analysis  concept,  for  large  numbers  of  modes,  N,  or 
large  differences  in  modal  frequency  (Aj  compared  to  XN),  time  domain  methods  are  numerically 
better  conditioned. 


2.  CHARACTERISTIC  SPACE  CONCEPTS 


A  new  way  of  conceptualizing  the  area  of  parameter  identification  was  developed  during  the  course 
of  performing  the  work  defined  under  this  contract.  One  of  the  objectives  of  this  contract  was  to 
summarize  existing  modal  parameter  estimation  methods  and  develop  new  ones.  In  the  process  of 
performing  this  task  it  became  obvious  that  most  of  the  current  algorithms  could  be  described 
conceptually  in  terms  of  a  three-dimensional  complex  space  of  the  system’s  characteristics.  Modal 
parameter  estimation  is  the  process  of  deconvolving  measurements  defined  by  this  space  into  the 
system’s  characteristics. 

The  frequency  and/or  unit  impulse  response  function  matrix  which  describes  a  system,  can  be 
expressed  in  terms  of  the  convolution  of  three  fundamental  characteristic  functions;  two  complex 
spatial,  and  one  complex  temporal.  The  spatial  characteristics  are  a  function  of  geometry  and  the 
temporal  corresponds  to  either  time  or  frequency.  Mathematically  the  frequency  response  matrix  and 
the  impulse  response  matrix  can  be  expressed  as  follows: 

[//(«*)]  =  [*]fA*J[L]  [A  </*)]-[*]  f '“mi  (22) 


where: 

•  [H (w*)!*,  x*,  *  frequency  response  matrix  (element  //„(w*)) 

•  I*  (* *)lw.  x*,  *  unit  impulse  response  matrix  (element  hM(tk)) 

•  I  *L».  xw  K  modal  vector  matrix  (function  of  spatial  variable  p,  element  j>w) 

•  [L  law  xtr,  *  modal  participation  factor  (function  of  spatial  variable  q,  element  L^) 

•  f  A*Jaw  xa*r  *  diagonal  martix  of  characteristic  roots  (element  - — - - ) 

I  JVk-K 

•  f  «WMa*r  xw  ■  diagonal  martix  of  characteristic  roots  (element 

•  wk*  frequency  temporal  variable  (k  *  1  — » biocksize  1 2,  may  be  unequally  spaced) 

•  tk  *  time  temporal  variable  (ik  *  k  A  t ) 

•  p  *  response  degree -of-frcedom  spatial  variable 

•  q  *  reference  degree-of-frecdom  spatial  variable 

•  r  *  temporal  degree -of-freedom  variable 

•  N  -  number  of  modes  (system  poles,  indexed  by  r) 

•  N„  *  number  of  responses  (indexed  by  p) 

•  *  number  of  references  (indexed  by  q) 

The  frequency  response  function  matrix  consists  of  three-dimensional  complex  space,  which  for  a  real 
system  is  a  continuous  function  of  the  three  characteristic  variables  (p,q,u> ).  However,  in  terms  of 
measurements  the  functions  consist  of  sampled  data  where,  p,q  and  wk  are  sampled  characteristic 
variables.  In  other  words,  the  frequency  response  function  is  measured  at  discrete  input,  or  reference 

I  points  ( q ),  output  response  points  ip),  and  discrete  frequency  (uk),  or  time  points  (tk). 

A  summary  of  the  characteristic  vectors  are: 


•  The  response  characteristic  functions  consist  of  a  set  of  vectors  which  are  proportional  to  the 
eigenvectors  of  the  system.  The  eigenvectors  are  indexed  by  r  and  the  elements  of  the  vectors  are 
indexed  by  p. 
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•  The  reference  characteristic  functions  consist  of  a  set  of  vectors  which  are  proportional  to  the 
modal  participation  factors,  which  are  in  turn  proportional  to  the  system  eigenvectors  at  the 
reference  degrees-of-freedom.  The  modal  participation  vectors  are  indexed  by  r  and  the  elements 
of  the  vectors  are  indexed  by  q. 

•  The  temporal  characteristic  functions  consist  of  vectors  which  are  equivalent  to  sampled  single 
degree-of-freedom  frequency  response  functions,  or  unit  impulse  response  functions.  The  index 
on  the  vector  is  r  and  the  index  on  the  sampled  element  of  each  vector  is  uk,  or  t  k. 

The  variable  r  is  the  index  on  the  characteristic.  For  a  given  r  there  is  a  discrete  characteristic  space. 
The  summation,  or  superposition  with  respect  to  r  defines;  the  measured,  or  sampled  frequency 
response,  or  impulse  response  matrix,  or,  in  other  words,  the  three-dimensional  complex  space. 

This  concept  is  difficult  to  visualize,  since  the  matrix  is  represented  by  three-dimensional  complex 
characteristic  space.  The  easiest  method  is  to  describe  the  variation  along  lines  parallel  to  axes  of  the 
space.  Lines  parallel  to  the  temporal  axis  correspond  to  individual  frequency  response  functions,  or 
unit  impulse  response  functions.  These  frequency  response  functions  consist  of  a  summation  of  the 
temporal  characteristics,  weighted  by  the  two  spatial  characteristics,  which  define  the  other  two  axis 
of  the  characteristic  space. 

Lines  parallel  to  the  response  axis  correspond  to  forced  modes  of  vibration.  These  forced  modes 
consist  of  a  summation  of  the  system  eigenvectors  weighted  by  the  input  characteristic  and  the 
temporal  characteristic. 

Likewise,  lines  parallel  to  the  input,  or  reference  axis  consist  of  a  summation  of  the  system 
eigenvectors  weighted  by  the  response  characteristic  and  the  temporal  characteristic.  The  variation 
along  these  lines  are  referred  to  in  the  literature  as  the  modal  participation  factors. 

Modal  parameter  estimation  is  the  process  of  deconvolving  this  sampled  space  into  the  basic 
characteristic  functions  which  describe  the  space.  In  practice,  there  are  many  more  measured,  or 
sampled  points  in  the  space  than  there  are  elements  in  the  three  characteristic  vectors,  therefore,  the 
parameter  estimation  process  is  over  determined.  As  a  result,  one  of  the  important  steps  in  the 
process  has  been  the  reduction  of  the  data  to  match  the  number  of  unknowns  in  the  parameter 
identification  process.  This  data  reduction  has  historically  been  done  by  using  a  pseudo  inverse,  or  a 
principal  component  method,  with  least  squares  being  the  most  common  pseudo  inverse  method. 

The  early  single  degree-of-fredom  (SDOF)  and  multiple  degree-of-fredom  (MDOF)  modal 
parameter  estimation  methods  used  subsets  of  the  sampled  data  and  extracted  one  of  the 
characteristic  functions  at  a  time,  normally  the  temporal  characteristic.  For  example,  the  very  early 
methods  like  the  complex  exponential  were  used  to  fit  individual  frequency  response,  or  unit  impulse 
response  functions  for  the  temporal  characteristics  (eigenvalue)  and  the  residues  (convolution  of  the 
response  and  input  characteristics).  For  these  cases,  each  frequency  response  measurement  gave  a 
different  estimation  of  the  system  eigenvalues,  or  temporal  characteristics.  Since  the  measurements 
were  taken  one  function  at  a  time  some  of  this  variation  was  due  to  inconsistencies  in  the  data  base 
and  the  rest  of  the  variation  due  to  noise  and  distortion  errors. 

Later  methods  started  to  use  either,  least  square,  or  principal  component  methods  to  condense  the 
data  over  a  number  of  sampled  frequency  response  functions,  into  small  subsets  parallel  to  the 
temporal  axis  (for  example  the  Least  Square  Complex  Exponential  and/or  the  Polyreference  Time 
Domain  methods).  These  methods  then  give  global  estimates  of  the  eigenvalues,  or  temporal 
characteristic  functions.  The  Least  Squares  Complex  Exponential  parameter  estimation  algorithm 
reduced  the  information  to  a  single  function  parallel  to  the  temporal  axis  and  as  a  result,  only 
estimated  the  temporal  characteristic  in  a  global  sense.  The  Polyreference  Time  Domain  algorithm 
estimates  several  functions  parallel  to  the  temporal  axis  at  the  input,  or  reference  points.  As  a  result 


this  method  also  gives  global  estimates  of  the  input  characteristic  functions,  or  modal  participation 
factors. 

The  more  recent  methods  use  larger  subsets  of  the  sampled  data  and  utilize  simultaneous  data  from 
all  three  axis  resulting  in  global  estimates  of  all  three  characteristics.  In  order  to  use  these  global 
methods,  it  is  important  that  a  consistent  data  base  be  measured. 
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3.  THEORY  OF  MULTIPLE  REFERENCE  PARAMETER  ESTIMATION 
TECHNIQUES 


3.1  Polyreference  Tune  Domain 


3.1.1  Introduction 

In  the  last  decade  modal  analysis  has  become  important  in  the  design  process.  There  are  two  rea¬ 
sons  for  this.  First,  there  have  been  significant  developments  in  computer  hardware.  More 
powerful  computers  have  made  the  heavy  computational  load  involved  in  modal  analysis 
achievable  within  reasonable  time  periods.  Furthermore,  developments  in  electronics  have  resulted 
in  larger,  cheaper  and  more  reliable  data  acquisition  systems.  In  addition  to  developments  in 
the  hardware  sector,  there  have  been  significant  developments  in  software.  New  parameter  esti¬ 
mation  algorithms  have  been  developed  to  provide  more  accurate  results. 

The  subject  of  this  section  is  one  of  the  latest  of  the  parameter  estimation  algorithms,  namely  the 
Polyreference  Time  Domain  Technique^18),  as  developed  by  H.  Void  at  S.D.R.C.. 

The  Polyreference  Time  Domain  Technique  is  a  two  stage  parameter  estimation  technique.  In  the 
first  stage,  the  damped  natural  frequency  and  damping  values  are  extracted  from  time  data  using 
single  or  multiple  references.  The  modal  coefficients  or  residues  are  then  calculated  in  the  second 
stage.  Although  the  first  stage  utilizes  a  time  domain  technique,  calculation  of  the  residues  can  be 
done  in  the  time  domain  as  well  as  in  the  frequency  domain.  When  only  data  acquired  from  a  single 
reference  is  used  to  estimate  frequency,  damping  values  and  the  modal  coefficients,  the  Polyrefer¬ 
ence  Time  Domain  Technique  is  equal  to  the  Least  Square  Complex  Exponential  l16l  Therefore, 
one  can  use  the  Polyreference  Technique  with  a  single  reference  or  with  multi-references. 

3.1.2  General  Equations 

The  Polyreference  Time  Domain  Technique  is  used  to  extract  modal  parameters,  namely  the  modal 
damping,  frequency  and  residues,  from  free  decay  or  impulse  response  function  data.  Instead  of 
utilizing  a  frequency  response  function,  like  many  frequency  domain  algorithms,  the  PolyTeference 
Time  Domain  Technique  is  formulated  from  the  impulse  response  function,  h(t).  This  h(t) 
represents  the  response  of  the  system  due  to  a  Dirac  impulse.  The  impulse  response  function  may 
be  obtained  by  taking  the  inverse  Fourier  transform  of  the  frequency  response  function. 

h(t)  =  Fl  (WW) 

Where : 

•  h(t)  =  impulse  response  function 

•  Pl  =  inverse  Fourier  transform 

•  H(ui)  =  frequency  response  function 

The  impulse  response  function  represents  a  decaying  time  history  of  the  system  and  can  be 
expressed  as  a  function  of  the  damped  exponentials  of  the  system. 

MO-EA  (23) 

rol 


Where : 


p  »  response  point 
q  -  input  point  or  reference  point 
r  »  mode  number 

-  residue  for  mode  r  at  point  p  due  to  an  input  at  point  q 
A*„ r  *  complex  conjugate  o fA^ 

Ar  *  system  pole 

Mffr*  /  H 

ar  *  system  damping  for  mode  r  (rad/sec) 

H.  *  damped  natural  frequency  for  mode  r  (rad/sec) 

A*  »  complex  conjugate  of  A, 

N  ■  number  of  modes  in  the  frequency  bandwidth  of  interest 


forr-lA3_N  {  f 
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Since  sampled  data  is  used,  h(t)  is  not  a  continuous  function,  but  a  discrete  function. 
w 

rm I 

Where : 

•  t„~kAt 

•  At  -  sample  time  or  time  resolution 

•  k  » integer  value 

•  K  gK*  A» 

•  2,  =  eA,A* 

By  substitution,  Eq.  (24)  becomes  : 

fsl 

IN 

<4>- (25) 

rat 

Because  there  typically  can  be  more  than  one  reference,  each  input  will  generate  an  equation  like 
Eq.  (24).  Therefore,  for  every  response  point,  a  set  of  N,  equations,  of  the  form  shown  in  Eq.  (26), 
can  be  written. 
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*pjfi-i(^*) =  E  ^pWi-if e ' 

r=l 

W'*)  =  E^,r 

r=l 


Where :  N{  -  number  of  references  or  inputs 

In  basic  modal  analysis  theory  I58!,  a  relationship  has  been  established  between  the  residues  and  the 
eigenvectors : 

Aptr  =  Qr  r  V>«r  r=  1,2,....,2N 
Where : 

•  Qr  =  scaling  factor  for  the  r-th  mode 

•  Vy  ”  scaled  p-th  response  of  a  complex  modal  vector  for  mode  r 

•  Vy  “  scaled  q-th  response  of  a  complex  modal  vector  for  mode  r 

This  relationship  enables  us  to  express  all  residues  as  a  function  of  the  residues  of  one  particular 
reference. 

(27) 

Writing : 

— (28) 

Vv 

where  is  defined  as  the  modal  participation  factor,  sometimes  referred  to  as  the  weighting  or 
forcing  factor,  i  is  a  reference  location,  q  the  reference  or  input  location  used  as  reference,  and  r  the 
mode  number. 

Substituting  Eqs.  (27)  and  (28)  into  Eq.  (26)  results  into  Eq.  (29). 

M4>-  E'W 

rsl 

^pa(4)  =  E^Zlr^plr  e*r  k 

rssl 


^pWi-l(Z*)=  E^-l.lr^pir  e 
f=sl 

W  JL  . 

W'*)  =  E^lr/tplr 

r=l 


These  equations  have  been  written  as  a  function  of  the  first  reference.  This  is  arbitrary  and  done 
for  convenience.  The  following  application  of  theory  is  independent  of  the  reference  used  for  the 
modal  participation  factors. 

Using  the  concept  developed  in  Eq.  (25),  Eq.  (29)  can  be  written  in  matrix  notation  : 

*pa(4)  I  r  /  X 


Where : 


: 


z1  0  0 
0  Z]  0 


0  0  0 


.2j=*A,A‘ 


-^pl.3 

=  residue  vector  = 


^3,1,1 
^  8*1,1 

^3,1,3 

. .  £-a,i,aw 

^W|l,3 

•  •  ^Nyl,2N 

•  \J-Hl  ](JV,-l)x2ff  =  ... 

•  (I  }lx2tf  =  a  row  vector  for  which  all  entries  are  1 

•[k]w,x2ff=  V  \  =  modal  participation  matrix 

1  1 


fV(4)1 

Mk) 


=  [l]  |ZJ‘{-4p} 


To  simplify  the  notation,  the  following  convention  will  be  used  : 


N,xl 


Equation  (30)  can  now  be  written  in  a  more  concise  form. 

I  **}  =  [^]*ixM  \Z\2Nxw  y4fj 

(  Iw,xi  v  Jawxi 

Note  that  only  the  left  hand  side  of  Eq.  (31)  is  known,  based  upon  measured  data. 


(31) 


3.13  First  Stage  Solution:  Poles  and  the  Modal  Participation  Matrix 

In  this  section  a  solution  will  be  developed  for  the  poles  and  the  modal  participation  matrix.  It  will  be 
shown  that  the  poles  are  the  eigenvalues  of  a  matrix  polynomial  and  that  the  columns  of  the  modal 
participation  matrix  are  the  eigenvectors  of  this  matrix  polynomial.  The  matrix  polynomial  is  the 
characteristic  equation  in  the  case  of  multiple  references. 

A  matrix  polynomial  is  a  polynomial  which  has  matrices  as  coefficients  instead  of  scalars  l59l.  The 
matrix  coefficients  of  the  polynomial  can  be  derived  in  different  ways.  In  this  section,  the  Prony- 
method  is  used  to  obtain  these  matrix  coefficients.  In  Section  3.1.5  the  matrix  coefficients  are 
derived  from  a  relationship  between  the  Laplace  transform  and  the  Z-transform. 

Using  a  matrix  polynomial  instead  of  a  scalar  polynomial  enables  the  Polyreference  Time 
Domain  Technique  not  only  to  process  data  from  multiple  reference  locations  simultaneously  but 
also  makes  use  of  the  fact  that  multiple  reference  data  is  available.  The  advantage  of  this  is  that  the 
Polyreference  Time  Domain  Technique  is  able  to  detect  repeated  roots  equal  to  the  number  of 
references,  or,  to  the  dimension  of  the  matrix  coefficients.  A  proof  of  this  is  given  in  Section  3.1.6. 

It  should  be  mentioned  that  data  from  all  references  can  be  used  for  the  pole  calculation.  However, 
certain  references  can  be  selected  to  emphasize  certain  properties  of  the  system.  In  the  previous 
section,  /V,  was  defined  as  the  number  of  references.  In  this  section  this  definition  is  still  valid, 
however  N,  does  not  have  to  be  equal  to  the  number  of  references  used  during  the  test.  Instead,  yV, 
is  equal  to  the  number  of  selected  references  chosen  to  be  included  in  the  calculations. 

The  primary  concern  in  any  parameter  estimation  technique  is  the  order  of  the  model.  For  the 
Polyreference  Time  Domain  Technique,  this  problem  can  be  translated  as  the  choice  of  2N  unknown 
poles  and  2N  unknown  mode  shapes,  where  2N  is  the  assumed  number  of  degrees  of  freedom  of 
the  system.  In  the  first  stage  of  the  calculations,  a  solution  is  sought  only  for  2N  unknown  poles. 

As  mentioned  before,  the  unknown  poles  are  found  to  be  the  eigenvalues  of  the  matrix  polynomial. 
The  number  of  eigenvalues  in  a  matrix  polynomial  is  equal  to  the  order  of  the  polynomial  times 
the  dimension  of  the  matrix  coefficients.  As  defined  before,  the  dimension  of  the  matrix 
coefficient  is  equal  to  the  number  of  references,  while  the  number  of  unknown  poles  is  2N.  There¬ 
fore,  the  order  of  the  matrix  polynomial,  n,  has  to  be  : 


Where : 


•  n  =  order  of  the  matrix  polynomial 

•  2N  *  number  of  degrees  of  freedom 

•  =  number  of  references 

Due  to  the  discrete  nature  of  the  sampled  data,  it  should  be  noted  that  n  must  be  defined  as  the 

IN 

smallest  integer  that  will  satisfy  Eq.  (32).  In  the  case  n  >  —  there  will  be  nxN,  eigenvalues 
found.  Since  nxN,  >  2N,  there  will  be  some  computational  poles. 

A  matrix  polynomial  that  yields  the  poles  and  the  modal  participation  matrix  is  sought.  This  poly¬ 
nomial  will  be  shown  to  be  of  the  form: 

[a  (0)|  [L]  \Z\*  ♦  [a  (1^  [L]  [Z J-1  ♦  •  •  •  +  [a  (n^  [l]  |ZJ°  =  [o]  (33) 

Where :  [a  (01w,xjv,  =  i-th  matrix  coefficient 
In  complete  matrix  notation,  Eq.  (33)  can  be  written  as: 

[o]ff|X(»+i)ff,  [^](»+i)w,x2Az  =  [^]wixaw  (34) 

Where : 

•  [q]  -  [[*  («)W<*  (/.  -1)1 - - [a  (0)] ]  JV|X(»+i)ft 

[L]  (Zf 

•  [°i-  ; 

[L  ]  \Z\ * 

'  J(»+i)Ar,xw 

One  way  to  prove  the  validity  of  Eq.  (33)  is  to  start  from  the  matrix  l  G  ]  in  Eq.  (34).  This  matrix 
has  2N  columns  and  ( n  +1  )Nj  rows,  where  : 

(n  +\)xNi  =  (”+l  )*Ni  =  2N  +  Ni 

The  columns  are  2 N  vectors  in  a  2N  dimensional  space.  These  ZN  vectors  may  be  independent 
with  respect  to  each  other,  but  do  not  have  to  be.  In  case  that  the  2 N  vectors  are  independent,  they 
form  a  base  for  a  subspace  of  the  space,  of  dimension  2N  +/V,.  In  other  words,  these  IN  vectors  span 
a  IN-  dimensional  space  which  is  a  subspace  of  the  global  IN  *  Ni  dimensional  space.  Since  this  glo¬ 
bal  space  is  spanned  by  2 N  *  Nf  independent  vectors,  the  2 N  independent  vectors  of  matrix  [  G  ]  can 
be  taken,  expanded  with  N,  vectors,  independent  from  these  2 N  vectors,  such  that  a  complete  base 
is  defined  for  the  total  space.  It  is  possible  to  define  the  Ni  additional  vectors  in  such  a  way  that 
they  are  orthogonal  to  all  2 N  vectors.  In  other  words,  the  dot  product  of  the  N,  vectors  with  the 
IN  vectors  is  zero.  This  is  stated  by  Eq.  (34). 
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[c]*,x(*+ipf,  [^](»*1WxM  =  [o]jv,x2W 


Where : 

•  [Q  ]  *  the  transpose  of  the  Nt  vectors 

•  [  G  ]  «  the  2N  vectors  of  the  subspace 


In  the  case  that  the  vectors  do  not  form  an  independent  set  of  vectors,  they  do  not  span  a  2 N 
dimensional  space,  but  a  Y  dimensional  space,  where  Y  is  the  number  of  independent  vectors.  For 
this  case,  there  still  exists  a  set  of  N{  vectors  which  are  independent  and  orthogonal  with  respect  to 
these  Y  vectors  and  Eq.  (34)  is  still  valid.  In  other  words,  if  the  2N  vectors  are  independent  or  not, 
A/,-  vectors  can  always  be  found  such  that  Eq.  (34)  is  satisfied. 


Using  this  property  of  the  (G  ]  matrix,  an  equation  can  be  derived  to  determine  the  matrix 
coefficients  [a  (()].  The  method  used  is  called  Prony’s  algorithm  [26],  This  method  makes  use  of 
the  matrix  coefficients  [a  (i)]  and  Eq.  (31)  to  derive  a  set  of  equations.  Equation  (31)  is  repeated 
here  for  convenience. 


jv}  -  [L]  PJ*  f*,} 


A  first  equation  is  obtained  by  taking  Eq.  (31)  at  time  (in  other  words,  k  equal  to  0+/i„  ),  and 
multiplying  this  equation  by  the  matrix  coefficient  [a  («)].  The  next  equation  is  obtained  by 
shifting  the  time  by  one  time  interval  (so  that  k  is  equal  to  n„+ 1),  and  multiplying  the  equation 
by  the  matrix  coefficient  [a(n-l)].  This  procedure  is  repeated  another  n-1  times  so  that  the  last 
equation  is  equal  to  Eq.  1311  at  time  t„  in_ .  multiplied  by  the  coefficient  matrix  [a  (0)].  Summing 
this  set  of  equations  gives  : 


[*(»:)  j'woj-  [«(«}  [l] 


j'wj  =  [ac«-^[L] 


(35) 


{wj  =  [«(1^  [l] 

[«(0]j  [«(0^[zj 


The  index  n«,  in  Eq.  (35),  is  used  to  indicate  an  arbitrary  origin  of  the  time  history.  This  enables  Eq. 
(35)  to  be  written  at  different  starting  points  for  the  same  time  history. 


Note  that  the  right  hand  side  of  Eq.  (35)  is  Eq.  (33)  multiplied  by  fZJ"*{/lp}.  Therefore,  the 
right  hand  side  of  Eq.  (35)  equals  zero. 


(36) 
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Equation  (36)  has  the  coefficient  matrices  [a  (/)]  as  unknowns.  It  should  be  emphasized  that  Eq. 
(36)  is  valid  only  for  free  decay  functions  or  for  impulse  response  functions  because  Eq.  (36)  has 
been  derived  from  Eq.  (31),  which  has  the  same  limitations. 

Since  [a  (0)]  is  the  matrix  coefficient  of  the  highest  order  in  the  matrix  polynomial,  [a  (0)]  can  be 
chosen  as  the  identity  matrix  [/  ],  without  a  loss  of  generality  f60).  Therefore  Eq.  (36)  can  be  written 
as : 


(37) 


For  each  point  p,  a  set  of  equations  can  be  written  by  varying  rta  from  0  to  (X-l),  where  X  is  an 
arbitrary  large  number.  An  identical  set  of  equations  can  be  obtained  for  every  response  location 
p.  All  of  these  sets  of  equations  can  be  written  as  one  equation  in  matrix  notation,  namely, 
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[flj  [[T?][T3  ] . .  •  mu]  =  -  [[/??]  W] . . .  [/?*„]] 


(38) 


Where: 

•  n  =  order  of  the  polynomial 

•  N„  =  number  of  response  locations 


•  [qj  =  [{*»«} .  { W ....  .{/w-i}]Kxx 

.  [flj-[[a(l)]  ,[a  (2)) . [a(n)] ] 

t»|X«W| 


{V-r} 

•  -  - 

•  .  .  {^rimJC-s} 

(W  {^l} 

•  •  •  {hpJC-l) 

•  X  =  large  number 


Note  that  the  last  row  in  the  [TJ]  submatrix  is  the  initial  portion  of  the  impulse  response  function, 
associated  with  each  respective  input  location.  The  second  to  last  row  is  the  same  impulse  response 
function,  shifted  over  one  time  sample  (A  /),  and  so  on,  as  shown  in  Figure  2. 


The  notation  can  be  simplified  by  using  the  following  definitions  : 

[a]  =  [(«(!)].[«  (2)]...  [«(«)]]  Ar,xJV,» 

[T].[mHT3]...m.)]w 


Equation  (38)  can  be  written  in  a  more  concise  form. 

[s]  [t]  =  [r]  (39) 

There  are  Af,x  nA/,  unknowns  in  Eq.  (39).  The  number  of  equations  is  N,  (  equals  the  number  of  rows 
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Figure  2.  [T]  -  Matrix 


of  [B  ])x  ( NoX )  (equals  the  number  of  columns  of  [T  ]).  By  taking  X  such  that  N„X  is  equal  to 
nNit  there  are  exactly  as  many  equations  as  unknowns  I60!  Therefore,  there  exists  a  unique  solution 
for  the  coefficient  matrices,  as  long  as  the  equations  are  consistent. 

However,  it  is  advisable  that  X  be  determined  such  that  N„X  is  greater  than  nN,.  This  makes  Eq. 
(39)  an  overdetermined  system.  Therefore  the  least  squares  method  can  be  used  to  solve  for  the 
coefficient  matrices.  Using  a  least  squares  method  provides  an  advantage  in  that  it  reduces  ran¬ 
dom  errors  in  the  data  t61l. 

However,  an  X  which  is  too  large  has  several  disadvantages.  Since  the  computational  load  is  not 
linear  with  this  variable,  it  is  not  advantageous  to  have  X  very  large  from  a  computational 
viewpoint  Second,  the  values  in  the  time  history  decrease  for  large  values  of  X,  since  the 
impulse  response  function  is  a  decaying  function.  As  a  result,  the  signal  to  noise  decreases  with 
increasing  time.  Third,  the  truncation  error,  due  to  an  inverse  Fast  Fourier  Transform  (F.F.T.), 
increases  at  the  end  of  the  data  block  (Figure  3).  Note,  this  truncation  error  only  exists  for  F.F.T. 
data. 


The  least  squares  solution  of  Eq.  (39)  can  be  found  from  the  normal  equation: 

M<£[U]  [<>--£  [«.][<  (40, 


where  the  entities  [TJ  ]  [TJJ  ]*  and  [/?£  ]  [TJ  ]r  can  be  shown  to  be  asymptotically  equal  to  the 
lagged  auto-  and  cross-  correlation  matrices  t62i  The  unknown  coefficient  matrices  can  be 
solved  from  Eq.  (40)  by  using  a  simultaneous  equation  solution  technique  such  as  Gaussian 
Elimination. 


Once  the  coefficient  matrices  have  been  estimated,  the  natural  frequency  and  damping  values  can 
be  determined  by  Eq.  (33). 
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Figure  3.  Impulse  Response  Function 
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Both  sides  of  Eq.  (41)  are  postmultiplied  by  a  unit  vector  (2 jVx  1),  composed  of  all  zeros  except  for 
unity  at  the  position  corresponding  to  the  system  pole  to  be  extracted. 


Where  : 

.  zr  •  (ek*At)*4 

•  {Ld  }/ir|Xi  =  column ’d'  of  [l  j 

•  {0  }W|Xl  =  zero  vector 

Since  {Ld  }  is  a  non  zero  vector,  Eq.  (42)  can  only  be  satisfied  if : 

det  [  fa  (/j|  zY  ]  =  0  (43) 

i=0  1  J 

Therefore,  the  resonance  frequency  and  damping  values  are  obtained  by  solving  for  the  poles  of  the 
matrix  polynomial. 

IN 

As  mentioned  earlier  this  matrix  polynomial  has  n/V,  =  ( — xNi)  =  2N  poles,  or  sometimes 

Ni 

(tt  ♦  1)  x  Af,  =  2iV  +  /V,  poles. 

Ni 

One  method  of  solving  this  system  of  equations  is  to  formulate  Eq.  (42)  as  a  standard  eigenvalue 
problem.  Recalling  that  \a  (0)]  =  [/  ],  the  following  matrix  polynomial  equation  is  obtained  from  Eq. 
(43). 


WWrMn 


I 


(44) 


j[j  (l)j  2J'1  +  [a  (2]j  zT2  + . . .  +  ]a  (n  jjx  jz^j  =  jz^j  23 


The  roots  of  Eq.  (44)  can  be  found  by  using  a  companion  matrix  approach  l63).  This  approach  uses 
(n-1)  matrix  equations  expressing  identity  relationships.  This  method  converts  the  polynomial  into 
an  eigenvalue  problem  of  the  form : 

([E]*^[/])  X[X  }={0}  (45) 

Where : 

•  0  =  eigenvalue 

•  {X  }  =  eigenvector 

•  [£  ]  =  the  coefficient  matrix 

Before  building  the  companion  matrix,  the  following  vectors  are  defined  to  enable  (n  -1)  identity 
equations  to  be  written: 

{F0}  =o«2S{Ld} 

{Vi}  = 

{ V 2}  =  <*23  1} 


}=  oy23'1  {Ld)=zd{V^) 
where  is  a  proportionality  constant 

Using  Eq.  (44)  and  Eq.  (46)  the  following  matrix  equation  can  be  formulated  as  shown  in  Eq. 
(47). 


T  -‘T 


[0]  [01 


[a (n-1)]  - (a(n) 

[0]  [0] 


[/]  [0] 


[K.i } 

f(^-i }' 

(V»*) 

=  zd 

{^0} 

{V0} 

kN,x1 

Equation  (47)  is  the  companion  matrix  equation  for  the  matrix  polynomial  Eq.  (44).  It  should  be 
noted  that  the  poles  of  Eq.  (47)  are  the  same  as  those  for  Eq.  (44).  The  vector  {Ld  }  in  Eq.  (44)  is 
the  eigenvector  associated  with  the  pole  zd  and  also  the  modal  participation  factor.  This  vector 
{Ld  }  is  proportional  to  { V0  },  as  defined  in  Eq.  (46). 

By  using  the  appropriate  solution  algorithm  t64!,  the  eigenvalues  and  eigenvectors  of  Eq.  (47)  can 
be  found.  The  resonance  frequency  and  damping  can  be  calculated  by  substitution  of  the  following 
equations : 


ss 


Although  the  pole  calculation  for  the  Polyreference  Time  Domain  is  a  purely  time  domain  technique, 
the  residue  calculation  can  be  performed  in  the  time  domain  as  well  as  in  the  frequency  domain.  The 
residue  calculation  in  the  time  domain  will  be  first  discussed,  followed  by  a  brief  discussion  of  the  fre¬ 
quency  domain  calculation. 


3.1.4.1  Time  Domain  Residue  Calculation 


To  determine  the  residues,  in  the  time  domain,  and  the  resulting  modal  coefficients,  Eq.  (31)  is 
utilized. 

jv}  “  [kjjv.xaw  [ZJawxaw  (49) 

v  Jw,xi  l  Jjwxi 


Matrix  [ZJ  is  a  diagonal  matrix  with  the  poles  on  the  diagonal,  with  every  pole  appearing  with 
its  complex  conjugate.  Therefore,  matrix  fZJ  can  be  rewritten  as : 


(50) 

awxw 

Where : 


Since  the  columns  of  the  [L  ]  matrix  were  obtained  as  eigenvectors  of  an  eigenvalue  problem 
that  had  z,  as  its  eigenvalues,  these  vectors  also  appear  in  complex  conjugate  pairs.  Therefore,  the 
[L  1  matrix  can  be  written  as : 

[l]  =  [[£],[£*]]  (51) 

Where : 


Zi  0  ..0 

0  z3  . .  0 
.  .  ..  0 
.  .  ..  0 
.  .  ..  0 
0  0  ..  z„ 


Concentrating  on  the  matrix  product  between  [L  ]  and  [ZJ,  this  matrix  consists  of  two  subma¬ 
trices  which  are  complex  conjugates  of  one  another. 


M  W- &*«•■**!]  [ft  |P|* 


Where : 


•  k  =  time  point 

The  residues  also  appear  in  complex  conjugate  pairs,  corresponding  to  complex  conjugate  poles. 
So  vector  [Ap  }  in  Eq.  (49)  can  be  rewritten  as : 


_  W 


Where : 


M,ul 

*  -^pia 


•H-M-  : 

ApiN 

>  4 

Substituting  Eqs.  (52)  and  (53)  into  Eq.  (49)  results  in  : 

r  i  [M 


R -2R'My’}} 


Where  :  Re  =  real  part  of  a  complex  number 


Note  that  V  and  U  are  complex  vectors. 


I/,  j  ■  jtV-wJ  *  i 


Where  : 


-  [K«,]-  Re|[K]| 

•  [^]=  Im|[K]J 


•  Re=  real  part  of  a  complex  number 

•  Im  =  imaginary  part  of  a  complex  number 


Equation  (54)  can  be  written  as : 


m  ■  2  [iv^  i* .  -  j*i  x{{(^)} 


In  certain  cases,  it  is  correct  to  assume  that  the  modes  are  normal.  In  these  cases,  {UPtnai }  ={0} 
l58!,  and  Eq.  (55)  can  be  simplified  as : 


A  least  squares  estimate  of  the  residues  can  also  be  formulated.  By  varying  k  in  Eq.  (49),  the  fol¬ 
lowing  equation  is  obtained: 


({V, }) 

[kpl  } 


[l  ]  rzj1 


(V)  \[L]\Z\X 

<  l  1  Jpr+i)JV|X2tf 


Where  :  X  =  an  arbitrary  positive  integer 


i  ‘  t  *  I.'  ' 


1 


Premultiplying  both  sides  by 


2*X 


[L‘  ]  r  it 
[L‘]  Z*  1 


results  in  the  normal  equation  as  follows : 


mH- mW 


Where : 


j^-jawxaw  =  [®]  x  ] 


{ftP t) 


'(X+iJjv.xi 


K  V ) 


From  Eq.  (58),  the  unknown  residues  can  be  determined  by  using  a  simultaneous  linear  equa¬ 
tion  solution  technique. 

A  least  squares  solution  can  also  be  found  for  the  normal  mode  assumption.  By  varying  k  in  Eq. 
(56)  a  set  of  equations  is  obtained,  similar  to  Eq.  (57) : 


J(*+i)w,  XI 


Where : 


pr+i)ffjxw 


:m{4 

[  t  Vwnag  ]°1 


By  premultiplying  both  sides  of  the  equation  by  the  transpose  of  the  coefficient  matrix  [D  ],  the 
normal  equation  is  obtained,  which  can  be  solved  again  by  a  simultaneous  linear  equation  solution 
technique. 


^wre "  *  [^]T  [^] 


3.1.4 2  Frequency  Domain  Residue  Calculation 


As  mentioned  before,  once  the  poles  and  modal  participation  matrix  are  known,  the  residues  do  not 
have  to  be  calculated  in  the  time  domain.  The  advantage  of  calculating  the  residues  in  the  frequency 
domain  is  that  the  effects  of  modes  outside  the  frequency  range  of  interest  can  be  compensated  by 
the  use  of  residuals.  The  disadvantage  of  the  frequency  domain  calculation  is  that  for  a  wide  fre¬ 
quency  range  (order  of  2  difference  between  the  lowest  and  highest  frequency),  this  calculation  is 
numerically  not  as  stable  as  the  time  domain  calculation. 


In  the  basic  modal  analysis  theory  l5®)  has  been  proven  that  the  frequency  response  function  in  the 
frequency  domain,  can  be  written  as : 


+ A, 


~)  +  Rr~  Rj 


where : 

•  Rr  *  residual  flexibility 


•  Rt  *  residual  inertia 


-  constant 


for  the  case  of  a  single  reference.  For  multiple  references,  the  previous  equation  becomes : 

|w>j»[L]  fC(w)J  |a,|  j  *  j 
where : 

•  *  modal  participation  matrix 


(61) 


•  iRr  }w,xi  *  flexibility  residue  vector 


•  iRt  }«,xi  =  inertia  restraint  vector 


lilLl 

u? 


•  fC(w)Jawxatf  * 


CnJI 


CN+l,N*l 


Where 

•  Cw,  =  •— —  for  i  =  1,2,.  JV 

jw\ 

•  =>  — l—  for  i  =  N +\,N +2,...,2N 

;w-A, 


-33- 


•  {^P}W(K  1  =  ' 

Equation  (61)  can  be  simplified  to  the  following  Eq.: 

-  [D]  f£(o»)J  j/?J  =  [>(«}  j*  J 

Where : 

•  [D  Uxapf-w,,  =  [[£]  \I\  f/j]  Where  f/J  is  the  identity  matrix 

rCJ  [0]  [0] 

•  r^Ja(#+*i)X3(Jir+j»i) =  I®]  f/J  [0] 

[°]  [0]  -j 

•  [FMbt&jfjf+fr,)  =  [D  ]  f£(w)J 


(62) 


•  iR  )HH*K i)xl  * 


Equation  (62)  is  a  set  of  Nt  equations  in  2 (N  +  A/,- )  unknowns.  In  order  to  have  a  unique  solution  for 
these  2(N  +  )  unknowns,  at  least  2(N  +  jV,- )  equations  are  needed.  If  more  than  2 (TV  +  N, )  equa¬ 

tions  are  available,  an  overdetermined  system  is  obtained  and  a  least  squares  estimation  can  be  used  . 

By  using  different  values  of  u  in  Eq.  (62)  an  overdetermined  system  is  obtained  of  the  form  : 


where  X  is  an  arbitrary  large  integer  which  has  to  be  smaller  than  the  used  blocksize  during  the  meas¬ 
urements. 


This  set  of  simultaneous  equations  can  be  solved,  for  example,  by  a  Gauss  elimination  process. 


It  should  be  noted  that  the  solution  for  Eq.  (64)  gives  complex  residues  while  the  residuals  are  real 
numbers. 

For  the  case  where  real  modes  are  desired  the  same  procedure  is  followed.  In  Eq.  (61),  the  matrices 
[L  ]  and  [C  ]  are  replaced  by  [Z/  ]  and  [C'  ]  where  : 


•  [C  Ww  = 


«u 

0  . 

.  0 

0 

c2,2  • 

• 

6 

6  ; 

•  CN,N 

Where : 


•  Ci.  = 


ju>  -  A, 


1 


V'w-A* 


•  [L  Iwixiir  =  Vi  l?  •  •  •  M 
[^Wm  1  [/1  I7  •  •  •  Ih  h  I7  •  •  ■ 

Performing  the  same  substitutions,  a  similar  equation  to  Eq.  (64)  is  obtained.  Solving  this  equation 
gives  real  modes. 


3.1.5  Time,  Frequency,  Laplace  and  Z-Domain  Relations 

In  this  section  the  validity  of  Eq.  (36)  will  be  proven.  The  proof  is  based  on  the  relationship 
between  the  Laplace  domain  and  Z-domain.  The  main  difference  between  these  domains  is  that 
a  Laplace  transformation  is  used  for  continuous  data,  while  the  Z-transformation  is  used  on 
sampled  data.  The  Laplace  Transformation  transforms  a  differential  equation,  in  the  time 
domain,  into  a  polynomial  equation  as  a  function  of  s  in  the  Laplace  domain.  The  Z-transformation 
converts  a  difference  equation,  in  the  time  domain,  into  a  polynomial  equation  as  a  function  of  z  in 
the  Z  domain. 

The  differential  equation  for  a  one  degree  of  freedom  system  with  single  input  is  : 

mx(t)*cx(t)*kx(t)=f(l)  (65) 

The  Laplace  Transformation  of  Eq.  (65),  assuming  zero  initial  conditions,  is : 

(ms2 +  cs+k)X(s)  =  F(s) 
or 


*SSlm _ l _ 

F(s)  ms2  +  cs+k 
Where : 

•  ms2 +  cs+k  =  Othe  characteristic  equation 

•  X(s)  =  Laplace  transform  of  the  response  x(t) 


(66) 


wwv.iV.'V-'/ x+i  vj  irjw w  wr j  *m  ^  rv  ^»r, « 


•  F(s)  -  Laplace  transform  of  the  force  /  (/) 

A  system  with  N  degrees  of  freedom  and  single  input  has  a  characteristic  equation  equal  to : 

^2S  +  fcjKj  ♦  bw-7  +  •  . .  +  bi  s  +  &o  =  0  (67) 

The  roots  of  this  polynomial  are  called  eigenvalues.  The  real  part  is  the  damping  and  the  ima¬ 
ginary  part  is  the  eigenfrequency  or  damped  natural  frequency. 

In  the  Z  domain,  z  is  equal  to  e' A*  with  A  t  the  sample  time.  Therefore  there  is  a  relation 
between  s,  the  Laplace  domain  variable  and  z,  the  Z  domain  variable.  Based  on  this  fact  it  can  be 
shown  that  there  exists  a  polynomial  in  the  Z  domain  that  has  exactly  the  same  roots  as  the  charac¬ 
teristic  equation  in  the  Laplace  domain,  namely: 


where  Ar  is  a  root  of  Eq.  (67).  This  new  polynomial  with  different  coefficients,  has  the  following 
form: 


+  02/r-i  z2*'1  +  ...  +a1z+a0  =  0 


(68) 


Assuming  now  that  the  different  characteristic  values  for  Eq.  (68)  are  known  and  can  be  represented 
by  Zr,  where  r  refers  to  the  root  number.  If  z,  is  substituted  into  Eq.  (68),  the  equation  is  satisfied  and 

can  be  multiplied  by  any  arbitrary  function  without  violating  the  equation.  Taking - —  as  an 

z  -  e  ' K 


arbitrary  function,  Eq.  (68)  becomes : 


(  «aw  2?  +  aaw-i  2™'1  *...+a12r*  a0) - ^r-r-  =  0  (69) 

2  -  e ’  * 

Where: 

•  tk=k&t 


Equation  (69)  represents  a  sum  of  constant  coefficients  multiplied  by  a  function  in  z.  Therefore,  if  an 
inverse  Z-transform  is  performed  on  Eq.  (69)  the  only  factor  that  changes  is  the  function  in  z  and 
Eq.  (69)  becomes : 

(flaw  *  flaw-1  z™’1  *  ■  ■  •  +  «1  2r  +  a0)  =  0  (70) 


or 

(  Efl,4)*Ar*  =  o 

•=0 

Substituting z,  =  ex’  A‘  into  Eq.  (71)  gives: 

£a<eM%*’A,)  =  0 
1=0 


(71) 


(72) 


Instead  of  — —  — ,  the  function  — PT  .  could  have  been  taken  as  arbitrary  function,  without  loss 
z-e's  2-erK 

of  generality.  For  this  arbitrary  function  Eq.  (72)  becomes : 


•'■Hi* 
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Equation  (73)  is  valid  for  all  characteristic  values  of  Eq.  (67).  In  other  words,  each  damped  exponen¬ 
tial,  associated  with  each  particular  characteristic  values,  satisfied  Eq.  (73).  Since  Eq.  (73)  is  a  linear 
equation  in  e,  any  linear  combination  of  the  individual  solutions  will  be  a  solution  for  Eq.  (73).  Recal¬ 
ling  now  the  equation  for  the  impulse  response  function,  namely : 

M'*)=  T,APreK%k 

fa O 

This  is  a  linear  combination  of  the  individual  solutions  of  the  linear  Eq.  (73),  and  is  therefore  also  a 
solution.  This  means : 

£**„(/*♦  f  A/ )-0  (74) 

•=o 

Where: 

•  A  t  =  sample  time 

•  h(fk)  =  impulse  response  function 

What  has  been  proven  up  till  now  in  this  section  can  be  summarized  as  follows.  Equation  (65) 
represents  the  differential  equation  for  a  continuous  function.  By  performing  a  Laplace  transform, 
this  differential  equation  is  converted  into  a  polynomial.  Then  by  showing  the  relation  between  the 
Laplace  domain  variable  s  and  the  Z  domain  variable  z  one  can  prove  that  a  polynomial  in  the  Z 
domain  exists  with  the  same  poles.  This  polynomial  is  the  Z  transform  of  a  difference  equation.  This 
difference  equation  is  the  equivalent  of  the  differential  equation  in  Eq.  (65).  Namely  a  differential 
equation,  which  is  only  valid  for  continuous  data,  converts  into  a  difference  equation  when  the  func¬ 
tion  is  a  discrete  function. 

It  is  known  that  the  solution  for  the  characteristic  equation  are  damped  exponentials  in  the  time 
domain.  For  the  case  when  discrete  data  is  available,  the  characteristic  equation  is  a  difference  equa¬ 
tion  in  the  time  domain.  Every  part  of  the  damped  exponential  satisfies  the  characteristic  difference 
equation.  This  property  has  been  proven  between  Eq.  (68)  and  Eq.  (73).  Namely  by  multiplying  by  the 

function  — z.  gives  in  the  characteristic  difference  equation  the  term  e*cV  By  taking  different 

z  -e  r% 

values  for  t,  different  portion  of  the  damped  exponential  function  are  used  and  the  characteristic 
difference  equation  is  still  satisfied.  This  proofs  that  indeed  any  portion  of  the  damped  exponential 
function  satisfies  the  characteristic  difference  equation  as  was  stated  previously. 

In  the  case  of  multiple  references,  Eq.  (67)  becomes  a  matrix  polynomial  of  the  form  : 

[/>  («)j  sn  +  [h(/i  -l]j  ♦  [f>  (rt  -2)|  s-3  ♦ .  . .  +  [b  (lj|  s  +  [ft  (0)|  =  [o]  (75) 

The  dimensions  of  the  matrix  coefficients  and  the  value  for  n,  the  order  of  the  matrix  polyno¬ 
mial,  have  been  discussed  in  Section  3.2. 

From  Eq.  (75)  and  using  a  similar  procedure  as  used  in  a  single  input  case,  the  following  equation  is 
obtained. 


iwtKBW  wognoror W  vtum*: 


or 
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nd 
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h(tk+i  At 


(76) 


D 


Where :  tk  =  k  At 

Restricting  k>  0,  Eq.  (36)  and  Eq.  (76)  are  identical. 

3.1.6  Detecting  Repeated  Roots 


The  single  reference  algorithms,  developed  prior  to  the  Polyreference  Time  Domain  Technique,  are 
theoretically  not  able  to  detect  repeated  roots  in  a  data  set.  However,  in  practice,  perfect  multiple 
roots  hardly  ever  occurs.  This  does  not  mean  that  very  close  coupled  poles  are  not  possible  in  a  data 
set.  Even  for  close  coupled  poles,the  single  reference  algorithms  may  have  trouble  to  detect  all  close 
coupled  poles  or  calculating  a  proper  estimate  of  the  modal  vectors.  When  a  modal  model  is  wanted 
for  later  applications,  such  as  modifacation  analysis  for  example,  missing  a  pole  or  having  a  poor 
modal  vector  estimate  in  the  modal  model,  makes  the  results  of  the  further  analysis  very  doubtful. 
The  Polyreference  Time  Domain  Technique  has  the  ability  to  detect  repeated  roots. 

This  proof  starts  with  Eq.  (36) : 

smM'H  <77> 

where  n.  is  taken  zero  for  convenience. 

As  mentioned  in  Section  2,  this  equation  is  only  valid  for  an  impulse  response  function  or  for  a  free 
decay  function. 


Consider  a  free  decay  function  of  only  one  frequency,  namely : 


Where : 


AP,N\,r] 


•  -  residue  for  mode  r  at  point  p  due  to  an  input  at  point  q 

•  A,  =  system  pole  =  ar  +jur 

•  ar  =  system  damping  for  mode  r  (rad/sec) 

•  =  damped  natural  frequency  for  mode  r  (rad/sec) 


(78) 


Since  Eq.  (78)  satisfies  the  conditions  for  which  Eq.  (77)  is  true,  Eq.  (77),  for  this  special  case, 
becomes : 
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Assume  now  that  the  pole  A,  is  a  repeated  root  with  multiplicity  two.  Further,  let  { A }  and  {Apr+ 1  } 
be  two  distinctive  residue  vectors  that  are  related  to  the  two  respective  eigenvectors  that  are  con¬ 
nected  with  this  pole. 

For  both  residue  vectors,  Eq.  (79)  becomes : 


t* 


Taking : 


(E[«(^xN,^,^)  =  [C(r)]WlXWl 


(80) 


Equation  (80)  can  be  written  as 


J^(r)]jv,  xn,  |j/^*rj  j/^*>'r+1 1  =  [^]v. 


X2 


(81) 


Therefore,  the  matrix  [C  ]  transforms  the  {Apr  )  residue  vector  as  well  as  the  {y4pr+1  }  residue  vec¬ 
tor  into  the  zero  vector.  As  a  result,  the  nullspace  of  the  transformation  matrix  is  at  least  two, 
since  {A^r  }  and  {A^x  }  are  independent. 

If  the  nullspace  has  a  maximum  dimension  of  at  least  two,  the  rowspace  has  a  dimension  of  max¬ 
imum  Ni  -2  with  Ni  the  number  of  rows  in  the  [C  ]  matrix  (Figure  4). 

Therefore,  to  determine  if  a  pole  has  a  multiplicity  greater  than  one,  the  matrix  [C  ],  for  this 
particular  pole,  can  be  calculated.  The  Tank  of  this  matrix  [C  ]  is  equal  to  the  dimension  of  the 
row  space.  The  dimension  of  the  nullspace  is  smaller  than  or  equal  to  the  multiplicity  of  the  pole. 
If  the  rank  of  the  [C  ]  matrix  is  different  from  zero,  then  the  dimension  of  the  nullspace  is  equal 
to  the  multiplicity  of  the  pole.  In  the  case  in  which  the  rank  of  the  [C  ]  matrix  is  equal  to  zero,  the 
multiplicity  of  the  pole  is  at  least  L,  the  number  of  rows  of  the  [C  ]  matrix.  Only  L  eigenvectors 
connected  with  that  pole  can  be  detected,  although  there  may  be  more. 
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3.2  Polyreference  Frequency  Domain 


3.2.1  Introduction 

Recent  advances  in  frequency  response  function  (FRF)  estimation  techniques  clearly  show  the 
benefits  of  multiple  input  I50-52).  The  redundancy  in  the  measurement  data  allows  the  analyst  to 
separate  closely  spaced  modal  frequencies,  to  identity  repeated  modal  frequencies  (repeated  roots)  or 
to  distinguish  between  closely  coupled  modes.  Noise  can  also  effectively  be  averaged  out  by  taking 
into  account  more  data  simultaneously. 

In  the  last  decade  several  new  techniques  have  been  developed  to  take  advantage  of  multiple  input 
measurements  such  as  the  Polyreference  Time  Domain  Technique  !18!and  the  Modified  Ibrahim 
Technique!65!.  These  methods  all  apply  to  time  domain  impulse  response  functions.  Since  the  actual 
measurements  are  commonly  stored  as  frequency  response  functions,  the  measurements  have  to  be 
transformed  to  the  time  domain  by  an  inverse  Fourier  transform.  This  inversion  gives  some  disad¬ 
vantages  to  the  time  domain  techniques: 

•  Truncation  errors  are  introduced  when  converting  data  from  frequency  to  the  time  domain  by  an 
inverse  Fourier  transform. 

•  In  order  to  be  able  to  perform  a  inverse  fast  Fourier  transform  (FFT),  the  data  has  to  be  taken 
with  a  constant  frequency  spacing. 

In  recent  years,  data  acquisition  systems  have  evolved  from  4  channel  systems  to  64  channels  or  more. 
This  development  minimize  the  disadvantage  of  sine  testing,  specifically,  the  time  consuming  nature 
of  the  sine  test  method.  This  is  due  to  the  fact  that  64  measurements,  or  more,  can  be  taken  simul¬ 
taneously.  Since  the  sine  technique  may  yield  superior  measurements  compared  to  random  excita¬ 
tion,  due  to  the  signal  to  noise  ratios,  there  is  a  tendency  to  evaluate  sine  testing  methods  again.  One 
of  the  advantages  of  sine  testing  is  that  the  frequency  spacing  does  not  have  to  be  constant.  How¬ 
ever,  existing  time  domain  modal  parameter  estimation  techniques  require  equidistant  frequency 
spacing  in  the  data  in  order  to  use  an  inverse  FFT.  To  take  optimal  advantage  of  sine  testing, 
without  resorting  to  a  forced  normal  mode  approach,  a  polyreference  frequency  domain  method  was 
necessary. 

With  this  background  in  mind,  a  Polyreference  Frequency  Domain  parameter  estimation  technique 
has  been  developed  simultaneously  at  the  University  of  Cincinnati  and  the  Katholic  University  of 
Leuven  (Belgium)  [42,43]  At  the  University  of  Leuven,  a  second  order  system  was  used  as  a  starting 
point  for  the  algorithm,  while  at  the  University  of  Cincinnati  a  first  order  system  was  used.  This 
makes  the  algorithm  developed  at  the  University  of  Cincinnati  basically  a  special  case  of  the  algo¬ 
rithm  developed  at  the  University  of  Leuven. 

In  this  section,  both  Polyreference  Frequency  Domain  Methods  will  be  described.  Both  methods  cal¬ 
culate  global  estimates  for  the  system  poles  as  well  as  the  modal  vectors. 

3.2.2  General  Concepts 

When  the  force  is  equal  to  a  unity  Dirac  impulse,  the  time  response  is  equal  to  the  impulse  response 
function.  This  means  that  this  time  function  can  be  expressed  as  a  sum  of  complex  exponentials  each 
multiplied  by  a  constant  (581.  This  constant  is  a  function  of  the  input  point  and  the  response  point. 
For  the  case  where  there  are  N,  inputs  and  N0  responses,  the  impulse  response  matrix  can  be  written 
as: 


5*? 
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[/«(/)]=  [*]e  W*[l]  (82) 

where: 

•  [A  (0  1/toxff,  =  impulse  response  matrix 

•  [♦  Ijvjxj#  =  matrix  with  the  modal  vectors 

•  fA  J wxw  ~  diagonal  matrix  with  the  system  poles 

•  [L  ]WxW|  =  modal  participation  matrix 

•  N0  =  number  of  measurements 

•  Ni  =  number  of  references 

•  N  =  number  of  system  poles 


B 


[L  ]  is  not  identical  to  the  modal  participation  matrix  as  it  is  defined  for  the  Polyreference  Time 
Domain.  The  scaling  factors  for  the  different  modal  vectors  are  included  in  this  matrix.  This  makes 
the  [L  ]  matrix  and  the  [¥]  matrix  together  a  unique  set  and  their  product  is  equal  to  the  residue 
matrix.  However,  the  same  interpretation  can  be  still  given  to  the  [L  ]  matrix  as  to  the  modal  partici¬ 
pation  matrix  as  it  is  defined  in  the  Polyreference  Time  Domain.  Namely  the  modal  participation 
matrix  indicates  how  well  a  particular  reference  excites  a  particular  mode.  Therefore,  the  [L  ]  matrix 
will  be  simply  referred  as  the  modal  participation  matrix. 

In  order  to  obtain  an  expression  for  the  displacement  transfer  function  in  the  frequency  domain,  the 
Laplace  transform  of  Eq.  (82)  can  be  taken,  which  gives: 

1,1  •  (83> 

where: 

•  s  =  Laplace  variable 

•  [H<t(s)  Ijvoxjv,  =  transfer  function  matrix  equal  to  the  Laplace  transform  of  the  impulse  response 
matrix  [h  (/)  ] 


•  Is  r/J  -  fAj  1  =  Laplace  transform  of  e 

L  J  2Nx2N 

•  f/J  =  identity  matrix 


:3 

V 

* 


The  subscript  d  on  the  [H  ]  matrix  indicates  that  the  transfer  function  matrix  is  in  terms  of  displace¬ 
ment  over  force.  By  defining  a  [  T  ]  matrix  as: 

[^jwxwi =  [*  r'J-  rAJ]  [^]  (84) 

Equation  (84)  can  be  simplified  to: 

[tf<r(s)]w^  =  [®]atoX2W  [tJwxw,  (85) 

In  order  to  obtain  an  expression  for  the  transfer  function  in  terms  of  velocity  and  acceleration  in  the 
frequency  domain,  the  single  and  double  partial  derivative  with  respect  to  time  can  be  taken  of  Eq. 
(82)  and  a  Laplace  transform  can  be  performed  on  the  obtained  equations.  Another  way  to  obtain 
these  two  equations  is  by  making  use  of  the  Laplace  transformation  properties  and  apply  these 
immediately  to  Eq.  (85).  The  resulting  two  equations,  obtained  for  the  transfer  function  matrix  with 
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respect  to  the  velocity  and  the  acceleration,  are: 


(86) 

=  s  [//„(*)]-  [•]  [l] 

(87) 

=  [*]  fAj  [r(s)j 

(88) 

and 

[*.<«>]  =  *3  [H*W]  •*  [h  <')«>]  -  [*(')«>] 

(89) 

=  s3[/fd(s)]-s  [*]  [l]-  f®]  fAJ  [l] 

(90) 

-  [*]  rAJ3  [ns)] 

(91) 

In  order  to  compensate  for  the  influence  of  the  modes  outside  the  frequency  range  of  interest,  an 
additional  term  can  be  added  to  the  previous  equations  to  compensate  for  these  residuals  effects.  Eq. 
(85),  (88)  and  (91)  become: 

[//«(*)]  -  [*]  [7(a)]  +  [/?„] 

[//.(*)]  =  [*]  taj  [m]  ♦  [r.] 

[«.(*)] » [•]  taj3  [m]  ♦  [ro] 

where: 

•  [ft*!  -  residual  effects  for  the  displacement 

l  Jwoxir, 

•  [/?,]  =  residual  effects  for  the  velocity 

l  iHoxN, 

•  Jr.]  =  residual  effects  for  the  acceleration 
l  JWoXWi 

3.2.3  Second  order  model 

In  this  section  the  more  general  second  order  model,  that  was  developed  at  the  Katholic  University 
of  Leuven,  will  be  discussed.  In  the  next  section,  the  special  case  of  this  model,  namely  the  first  order 
model,  will  be  discussed. 


(92) 

(93) 

(94) 


Equations  (92),  (93)  and  (94)  can  be  combined  into  one  equation: 


(95) 


Hd(s) 

H,(s) 

H.(s) 


M 


iWoxW 


SNoXN, 


Substituting  Eqs.  (87)  and  (90)  into  Eq.  (95)  gives: 

[Hd(s)]  m 

s  [//<(*)]  *  *]  [L] 

5a  [has)]-*  [*][l  -  [«]  fAJ  [l] 


[*1  M 


Concentrating  on  the  [[’*'],  (’*'  ]  fA  J ,  [¥J  [A  J2  ]r  matrix,  this  matrix  has  3N0  rows  while  2N 
columns.  This  means  that  this  matrix  represents  2N  vectors  in  a  3 N0  dimensional  space.  When 
3N0  >  2N  these  vectors  define  a  subspace  of  a  3N0  dimensional  space.  This  means  that  at  least 
3 No  *  2 N  linear  independent  vectors  can  be  found  that  are  orthogonal  to  these  2 N  vectors.  If  the  2 N 
vectors  are  not  independent,  more  than  3N0  ■  2N  orthogonal  vectors  can  be  found.  If  the  indepen¬ 
dent  orthogonal  vectors  are  represented  by  a  matrix  [Q],  the  previous  property  can  be  expressed  in 
next  equation: 


(3tfo-2tf)xSWa 


W-H 


where: 


•l©]3WoX2ff=  I1*)  fA  i 

[*]  TAJ2 

Under  the  restriction  3 N0  -2N>  N0m  or  N0  >  N,  at  least  N0  orthogonal  vectors  can  be  found  and 
Eq.  (97)  becomes: 


[fl]Z*».[e]  =  [o]A 


The  reason  for  the  restriction  N0  >  N  is  that  a  matrix  polynomial  is  desired  that  will  have  poles  equal 
to  the  system  poles.  This  restriction  for  a  matrix  polynomial  is  necessary  so  that  the  coefficients  will 
be  square  matrices. 

The  fact  that  the  different  vectors  in  the  [0  ]  matrix  appear  in  complex  conjugates,  makes  the  orthog¬ 
onal  vectors  monophase  vectors.  A  monophase  vector  is  a  vector  that  can  be  scaled  such  that  all  ele¬ 
ments  of  the  vector  are  real  or  imaginary.  An  example  of  a  monophase  vector  is  the  normal  mode 
vector.  The  monophase  property  can  be  proven  as  follows: 

Assume  that  the  complex  vectors  of  the  [0  ]  are  written  as 


H  ■  N  *'  M 


and  their  complex  conjugate  vectors  as: 


M'  ■  M  N 


where  j  *  v^-T 


Writing  the  orthogonal  vector  as  the  sum  of  real  plus  imaginary  part: 


In  order  for  the  [Q  ]  matrix  to  be  orthogonal  to  the  complex  conjugate  vectors,  it  has  to  satisfy  next 
two  sets  of  equations: 

(  [f?r]  *j  [Q<]  )*  (  [^]  +j  [O;] )  =  [o] 

and 

([a]+;‘N)x(N'-/  W)=[°] 

These  two  sets  of  equations  can  be  reduced  to  the  following  constraints: 


Qr 

= 

0 

Qr 

. 

et 

= 

0 

Qi 

= 

0 

Qi 

P 

xs 

0 

In  order  to  satisfy  these  restrictions,  the  orthogonal  vectors  have  to  be  monophase.  Practically,  this 
means  that  these  orthogonal  vectors  might  be  complex,  but  a  scaling  factor  can  always  be  found  such 
that  this  scaling  factor  makes  them  real  vectors.  Mathematically,  this  comes  down  to: 

M-H 

where: 

•  [a]  is  a  real  matrix 

•  a  is  a  complex  scaling  factor 

The  matrix  [Q]  can  be  divided  into  submatrices: 


B0 


Without  affecting  the  generality  of  Eq.  (98)  and  taking  advantage  of  the  monophase  property,  the 
matrix  [Q]  can  be  normalized  with  respect  the  submatrix  [B2\  and  the  [Q]  matrix  becomes: 


(99) 


Substituting  Eq.  (99)  into  Eq.  (98)  gives: 


t  nt  hr  wi 


l^ox^o  |  jj^I  r[Aji 


laWoxaw 


Equation  (100)  represents  a  matrix  polynomial  of  order  2.  Since  the  matrix  coefficients  have  a 
dimension  of  N0  this  matrix  polynomial  has  2N0  poles  The  IN  system  poles  represented  by  the  matrix 
fAJ  are  a  subset  of  these  2N0  poles.  Since  N0  can  be  a  lot  larger  than  2 N,  this  can  cause  a  numerical 
problems.  This  numerical  problem  can  be  minimized  and  this  will  be  discussed  later  in  this  text. 

Bringing  the  residue  matrix  in  Eq.  (96)  to  the  left  hand  side,  the  next  equation  is  obtained: 


[//«*(*)]  -  [Rd] 
s  [//<(*)  -  [*]  Z,]-  [/?. 
s2  [//„(*)]-*  [*]  [L]-  [*]  fAJ  [L  -  [*„] 


[ns)]. 


3H„X.2N 


l3WoXAf, 


Postmultiplying  Eq.  (100)  with  the  [T  ]  matrix  and  substituting  Eq.  (101)  into  the  obtained  equation 
gives: 

[//„(*)]  -  [Rd] 

[  M  Kf  f/J]  4  [^)]  *  *]  L]‘  [*.]  =  [o]  (102) 

s2  [//„(*) ]  -s  [*-]  [L  -  [*]  fAJ  [l]  -  [Ra] 

This  equation  means  that  the  measured  FRF’s  can  be  described  by  a  linear  frequency  domain  model 
with  real  matrix  coefficients: 

h]T  h(4)]+s  hf  h(4)]+s2  h<4)] = 

[  hr  M +s  W +  M rAj]  N *  hf  M +  hf  M +  M  (io3) 

~s  h] +  hi 

where: 

•  Mi]  3  !♦]  [L]  and  is  equal  to  the  residue  matrix. 

•  Mo]  =  [Mif  («]  +  [#]  [Al]  [L]  +  f/io  i  T[Rd  J  +  [A !  ]T  [/?.  1  ♦  [/?.  ] 

This  equation  holds  for  each  value  of  s  ( s  =jw  for  frequency  response  function  )  of  the  frequency  band 
where  Eq.  (82)  is  valid.  It  defines  a  linear  model  for  the  measured  FRF’s  between  N0  response  loca¬ 
tions  and  Ni  inputs  or  references.  The  unknowns  of  this  model  are  the  matrix  coefficients  \An  1 
Mil.  Mi]  and  Mo]- 

To  summarize,  the  algorithm  proceeds  as  follows:  first  the  unknown  matrices  [/t0]  ,  f/lj],  [/?,]  and 
Mol  are  calculated  by  solving  Eq.  (103)  in  a  least  squares  sense,  then  the  solution  of  the  generalized 
eigenvalue  problem,  defined  by  Eq.  (100),  directly  yields  the  complex  system  poles  and  the 
corresponding  mode  shapes  for  the  structure. 


As  already  mentioned,  Eq.  (100)  is  a  matrix  polynomial  that  has  2 N0  poles.  In  the  case  where 
N0>>N,  the  solution  of  the  eigenvalue  problem  of  dimension  2 N0  will  generate  a  lot  of 


computational  modes,  since  only  2N  physical  modes  ( including  the  complex  conjugate  system  poles  ) 
contribute  to  the  measurement  data.  Moreover,  the  solution  process  for  [>40]  ,  [/IjJ,  [flj]  and  [/?0] 
may  then  become  ill  conditioned.  To  solve  this  problem,  the  data  can  first  be  condensed  using  princi¬ 
ple  component  analysis.  I66*67! 

3.2.4  First  order  model 

This  section  explains  the  first  order  model,  which  is  a  special  case  of  the  second  order  model  which 
was  discussed  in  previous  section.  The  first  order  model  was  developed  and  implemented  at  the 
University  of  Cincinnati. 


a 


Instead  of  combining  Eqs.  (92),  (93)  and  (94)  to  develop  a  second  order  model,  Eqs.  (92)  and  (93) 
can  be  combined  into  one  equation  to  develop  a  first  order  model: 


Substituting  Eq.  (87)  into  Eq.  (104)  and  bringing  the  residual  matrix  to  left  hand  side: 


(104) 


j 


[«<(*)]  -  [*<] 

*  [HAs)]  -  [*]  [l]  - 


IWoXjV, 


[7’(s)]awxWi 


(105) 


Through  the  same  reasoning  as  in  the  previous  section,  it  can  be  shown  that  if  N0  >  2N  there  are  at 
least  N0  vectors  that  are  orthogonal  with  respect  to  the  2 N  vectors  in  [l*] ,  [*]  [a]1t.  This  property 
can  be  expressed  in  Eq.  (106)  where  the  [B  ]  matrix  containing  the^  orthogonal  vectors  is  already 
expressed  in  terms  of  submatrices: 

[*o 

Bx 


Since  the  2 N  vectors  in  [  [♦  ]  ,  [  V  ]  [AJ  ]  appear  as  pairs  of  complex  conjugate  vectors,  the  orthogo¬ 
nal  vectors  are  monophase  vectors.  The  matrix  [B  ]  can  be  normalized  with  respect  to  [AJ  and  Eq. 
(106)  becomes: 


Equation  (107)  is  a  first  order  matrix  polynomial  where  the  matrix  coefficients  are  real  coefficients 
with  dimension  N0.  Therefore,  the  matrix  polynomial  has  N0  poles.  Notice  that  the  matrix  polynomial 
for  the  second  order  model  was  of  order  2  with  2N0  poles,  which  could  lead  to  numerical  unstability. 
Therefore,  by  chosing  a  first  order  model,  this  problem  is  minimized. 

Postmultiplying  Eq.  (107)  with  the  [T]  matrix  and  substituting  Eq.  (105),  the  obtained  equation 
becomes: 


■5 

•j 
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j 

j 

V 

V 

i 


-47- 


! 


mmuwiwww  n  m  *.■  vv-V-VuV.y.^jyyvyvYr^  )prjnirtjr>jnfjrty*j-*jr*jr*jr*jrm. 


Ihf 


*oX2tf0 


[«<(*)]  -  [*d] 

M-MH-N 


*<>X*| 


2*oX*i 


or 


Mr  h(,)] =*4  h<s>] ♦  M  H  ♦  [*•] ♦  hf  M 

- -i  h^)]  ♦  h] 


where:  [B,  ]  =  [*]  [L  ]  ♦  [/?,  ]  ♦  [/l0  ]r  [*d  ] 


(108) 


(109) 


This  means  that  the  measurement  data  [//]  can  be  described  by  a  model  which  is  linear  in  the  real 
unknowns  [X0]  and  [Bt  ]. 

It  is  possible  to  reduce  the  dimensions  of  the  matrices  by  taking  advantage  of  some  system  properties. 
It  is  known  that  elements  of  ],  fA  J  and  [L  ]  occur  in  complex  conjugate  pairs.  This  property  can 
be  used  advantageously  to  reduce  the  dimensions  of  the  different  matrices  used  in  the  first  order 
model. 


If  the  complex  conjugate  terms  are  ignored  in  the  matrices  [*  ],  [A  J ,[T  ]  and  [L  ],  the  dimensions 
reduce  from  2N  to  N  and  the  reduced  matrices  can  be  represented  by: 


[*'] 

l  1*0 x* 


KL»  trwL», 


Using  this  notation  Eq.  (105)  becomes: 


\m4(s)]  -  [r;] 

•  KIM  -  N 


2*0X*i 


[*  ]  , 
[*]  KJ 


(HO) 


2*0 X* 


where  the  matrices  [Rd  ]  and  ]  no  longer  only  compensate  for  the  residuals  effects  of  the  out-of- 
band  modes,  but  also  compensate  for  the  error  that  is  introduced  by  neglecting  the  complex  conju¬ 
gate  term  in  the  different  matrices. 


The  matrix  [  [¥  1  [♦  ]  |a  J  ]r  has  only  N  vectors  in  a  2N0  dimensional  space.  Therefore,  when 
No  >  N,  there  exits  No  vectors  that  are  orthogonal  to  these  N  vectors.  This  is  expressed  in  Eq.  (Ill): 


(111) 


The  orthogonal  vectors  can  be  normalized.  In  this  case,  Eq.  (Ill)  becomes: 


\M 

[  [*•]  1 

VI 

h  [*l] 

(112) 


The  disadvantage  of  this  reduction  in  the  matrix  size,  when  the  complex  conjugate  terms  are  ignored, 
is  that  the  matrix  [A0  ]  is  no  longer  a  real  matrix  but  a  complex  matrix.  Postmultiplying  Eq.  (112)  by 
the  [7 ']  matrix  and  substituting  Eq.  (110)  yields: 
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jT*,jkL  jft* Mu  Jrfb  jIL  j£* 


ikft/j  [  M-M 

Ihl  f,)U».|I[H4W].[,]jL].[R;] 

or 

plof  [H4(s)]  =  -s  [//„(*)]  *  [tf  ]  [l‘  ]  ♦  ^of  [ 

•*  [«,(f)]  ♦[«;] 

where: 

•  [«;]-[*'  nr  ]+[A'0  ]♦[*:] 

This  means  that  the  measurement  data  [H  ]  can  be  described  by  a  model  which  is  linear  in  the  com¬ 
plex  unknowns  [A'0  ]  and  [B\  ]. 

3.2.5  Principle  Component  Reduction 

In  order  to  reduce  the  size  of  the  matrices  involved  in  the  Polyreference  Frequency  Domain 
Methods,  the  measured  data  can  be  transformed  to  their  principle  components  by  means  of  a  linear 
transformation!40'66'67].  This  process  minimizes  not  only  the  computer  requirements  (storage 
memoiy  and  execution  time)  but  also  provides  a  good  way  for  estimating  the  number  of  modes  con¬ 
tributing  to  the  response  of  the  structure. 

The  measurement  data  [//]  and  the  principle  components,  represented  by  [//*].  are  interrelated  by 
following  pair  of  linear  transformations: 

=  |cjw<xw<>  [//(5)]w#xw,  (115) 

or 

["<*)],.**  =  [c]W  (116) 

where  N,  is  the  number  of  principle  components  and  is  in  general  much  smaller  than  the  number  of 
measurements  N0. 

The  principle  component  data  [//']  will  satisfy  the  same  linear  model  as  the  actual  measurement  data 
[//].  Indeed,  combining  Eq.  (82)  and  (115)  results  in: 

M’[c]MrW‘H  <1,7) 

=  [**]e  [L] 

Therefore,  [//*]  will  also  satisfy  the  second  order  linear  model  with  real  coefficients  as  stated  in  Eq. 
(103): 


-  [0]*ox*.  0») 

ViaxN, 

**]♦[*;]  (114) 


(118) 


v.vv.v. 


[*o]  [hc(s)]+s  j^]  [//'(*)]  [//e(s)]=S  [*f  ]♦[/?§] 


where: 

•  i  R\)=m  [ L ] 

•  [*S1  =  [[Ax]  [V]  ♦  m  rAj][L]  ♦  [A0  ]  [R'd  ]  ♦  [Ax  ]  [1?;  ]  ♦  [ R .  ] 

The  same  algorithm  can  be  applied  to  the  condensed  principle  component  data  to  calculate  the 
modal  parameters.  The  results  for  these  data  will  be  a  set  of  eigenvalues  fAJ  and  eigenvectors  [4*]. 
Note  that  from  Eq.  (82)  and  (100),  fAJ  remains  the  same  for  the  principle  component  data  The 
mode  shapes  [4  ]  can  be  calculated  from  the  mode  shapes  [4*  ]  of  the  principle  component  data  by 
the  following  transformation: 


[c]w  [**], 


3.2.6  The  Modal  Participation  Matrix 

The  calculation  of  the  modal  participation  factors  [L]  of  Eq.  (82)  can  be  performed  on  the  actual 
data  as  well  as  on  the  principle  component  data.  The  system  of  equations  to  be  solved  is: 


[W(j)]^w,  *  [5  f'J  *  rAj]WxW  [L]awxw, 


where  the  mode  shapes  [4  ]  and  the  poles  fA  J  are  already  known.  To  keep  the  modal  participation 
factors  global,  all  output  stations  ( or  all  principle  components)  should  be  included  in  the  calculation. 

However,  the  matrices  [fig]  and  [/ff],  for  the  case  when  a  principle  component  reduction  is  per¬ 
formed,  contain  sufficient  information  to  solve  for  the  modal  participation  factors,  since  the  matrices 
fA  J  and  [4  ]  are  known  in  this  phase  of  the  parameter  estimation.  From  Eq.  (118)  the  matrix  [L  ] 
can  be  expressed  as  a  function  of  the  known  variables: 


4*]  f 

[*x]  [**]  ♦  [f]  fAJ  Rh] 


To  investigate  the  accuracy  of  the  parameter  estimation  process,  Eq.  (84)  can  now  be  evaluated  in 
order  to  curve-fit  the  measurement  data  for  some  or  all  output  and/or  input  locations.  The  errors 
between  the  actually  measured  and  the  synthesized  frequency  response  functions  will  be  minimum  in 
a  global  linear  least  squares  sense. 
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3.3  Multiple-Reference  Ibrahim  Time  Domain 
33.1  Introduction 


1 


The  Ibrahim  Time  Domain  technique  was  originally  developed  for  the  analysis  of  single¬ 
input/multiple-output  free  decay  or  impulse  response  time  domain  data  I17-28!.  This  technique  can  be 
expanded  into  the  multiple-input/multiple-output  algorithm  by  introducing  the  relation  of  mode 
shapes  and  modal  participation  factors. 

Principal  Component  Analysis,  Minimal  Realization  Approach  and  the  generalized  inverse  technique 
have  been  found  to  be  very  useful  during  the  development  of  the  Polyreference  Frequency  Domain 
l42-43!  and  the  Eigensystem  Realization  Approach!34].  Principal  Component  Analysis  is  used  to  evalu¬ 
ate  the  rank  of  an  auto  correlation  matrix  and  to  aid  in  the  selection  of  the  number  of  modes. 
Minimal  Realization  Algorithm  is  used  to  reduce  the  size  of  the  system  matrix  of  an  eigenvalue  prob¬ 
lem  to  the  approximate  rank  of  the  matrix.  The  generalized  inverse  technique  is  used  to  solve  the 
equation  [B  ]{*}  =  {b}  with  or  without  the  normal  equations  method.  Since  these  mathematical 
approaches  are  also  used  in  the  development  of  multiple-reference  Ibrahim  Time  Domain,  this 
method  is  similar  to  the  Polyreference  Frequency  Domain  and  the  Eigensystem  Realization  Algo¬ 
rithm.  It  can  be  said  that  Multiple-Reference  Ibrahim  Time  Domain  is  a  time  domain  version  of 
Polyrefetence  Frequency  Domain  or  another  approach  of  Eigensystem  Realization  Algorithm. 
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The  mathematical  derivation  starts  from  the  time  domain  equations: 
[/.(/)]-[♦]  [L] 

where: 

•  ['i/  ]  -  by  2N  mode  shape  matrix 

•  [  L  ]  =  2N  by  Ni  modal  participation  factor  matrix 

•  N  -  number  of  modes 

•  N„  *  number  of  outputs 

•  Nf  =  number  of  inputs 


Expanding  these  matrices  with  respect  to  sampling  times: 


=  m  [L] 

where: 


•  [/4(0]  = 


r  ih{t)] 
[h(t+&t)] 


IMfAt)] 

[h(t+2At)  ] 


[A(r+(i-l)Af)]  [A(f+/Af)l 


[  h  (t  *(k  -1)A  t)  ] 
[n(t+k&t)  ] 

[  h  (t+(i  +k  -2)A  t)  ] 


(122) 


(123) 
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[*] 

[*] 

[*] 


.[L]*[[L]  I  r*ArA‘j [L ]  I  ...  I  r^(*-^‘j[L)] 

•  i  =  number  of  time  shifts  in  column  direction 

•  k  »  number  of  time  shifts  in  row  direction 

The  matrix  [  //*(f)  ]  is  called  Hankel’s  matrix.  Note  that  the  matrices  [  //*(/)  ]  and  [  ♦  ]  are  nor¬ 
mally  rectangular  matrices,  therefore,  only  their  generalized  inverse  exits. 


Using  one  sampling  time  shift  on  Eq.  (123): 

[H*(f+Af)]*[*l  r<ArA‘j  r«MJ  IM 

Premultiply  Eq.  (123)  by  [  ♦  ]*  ,  the  generalized  inverse  of  the  [  ]  matrix. 

[•riflwoi-  r«Mj  m 


(124) 


(125) 


The  right  hand  side  of  this  equation  is  a  common  term  in  Eq.  (124).  Substituting  Eq.  (125)  into  Eq. 
(124)  yields 


[//*(/ +Af)]  =  [*l  re*,A‘j 

Denoting  [  *  ]  r  eArA‘ j  [  *  J*  as  matrix  [  B  ] : 
[J/*(/+A/)]»[£  1[//*(01 


(126) 


(127) 
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This  is  a  matrix  homogeneous  difference  equation  which  also  can  be  derived  via  Prony  s  algorithm. 
The  matrix  [  B  ] ,  is  the  matrix  of  coefficient  matrices  of  the  matrix  polynomial  equation. 

To_obtain  the  eigenvalue  _problem,  postmultiply  Eq.  (126)  by  the  generalized  inverse  of 
[  [  *  1*  [  W*(t)  1 1 ,  since  l  [  *  ]*  [  H*(f ) )  ]  is  in  general  not  a  square  matrix. 

[//*(/ ♦A/)]  [[*N //*(/))]♦  =  [*]  feArA‘j  (128) 

Using  the  property: 

Equation  (128)  can  be  simplified  to  : 

[B  ](*]  =  [*]  r«ArA*j  =  [*]  r«*rj  (129) 

where: 

•  [*]-["•('+*')  ][//*(/)  Y 

Equation  (129)  is  an  eigenvalue  problem;  the  size  of  the  system  matrix  [B  ]  is  N0i  xN„i  .  The 
eigenvalues  that  are  found  are  of  the  form  zr  =  eArA‘. 
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Mode  shapes  are  calculated  as  eigenvectors  of  Eq.  (129).  The  eigenvalues  are  translated  to  damped  \ 

natural  frequencies  and  damping  factors  as  follows  (  A,  =  a,  +  u,  ):  j 

*  Re  ( In 2, )  /  At  u,  =  Im  ( In 2; )  /  At  (130)  J 


3.3.2. 1  Minimal  Realization  Approach 


The  minimal  realization  approach  is  a  technique  that  is  used  to  reduce  the  size  of  the  system  matrix 
[  B  ]  of  an  eigenvalue  problem  to  the  approximate  rank  of  the  matrix. 


The  matrix  [  //*(f )  ]  is  decomposed  as  follows: 

[ff>k(f)]HJxNlk  =  [P  ]w„ixW»i  N„ixN,k  [  Q  ]^,»  xM* 

Denoting  the  rank  of  the  matrix  [  H^ft)  ]  as  p  ,  the  matrix  ( D  j  is  written  as  follows: 


di  0  .  .  0  0  .  .  0 

0  d2  .  .  0  0  .  .  0 


fDj  = 


0  0  .  .  dv  0  .  .  0 
0  0  .  .  0  0  .  .  0 


0  0  .  .  0  0  .  .  0 


(131) 


(132) 


These  diagonal  terms  are  called  singular  values.  The  number  of  significant  singular  values  and,  there¬ 
fore,  the  rank  p  of  the  matrix,  may  be  evaluated  using  the  definition  of  condition  number  and/or  the 
ratio  of  successive  singular  values.  This  procedure  is  explained  as  follows: 

Condition  numbers  are  defined  as  follows: 


CN  = 


dx 

dk 


(133) 


Since  the  first  singular  value  is  the  largest  numerical  value,  the  condition  number  is  always  larger  than 
zero.  If  the  condition  number  is  evaluated  for  different  submatrices  of  the  [D]  matrix  by  taking  suc¬ 
cessive  singular  values,  the  rank  p  of  the  matrix  is  defined  when  the  condition  number  first  becomes 
very  large. 


The  rankp  can  also  be  determined  from  the  ratio  of  successive  singular  values  calculated  as  follows: 
dh 

ratio  =  - —  (134) 

ak* 1 

If  the  ratio  of  successive  singular  values  is  evaluated,  this  ratio  will  be  approximately  equal  to  one 
when  the  successive  singular  values  are  significant.  This  will  also  be  true  when  the  successive  values 
are  trivial.  At  the  transition  between  these  two  subsets  of  the  singular  values,  the  ratio  of  successive 
singular  values  will  peak  indicating  the  rank  p  (  p  =  k  ).  Then,  the  dimensionality  of  Eq.(131)  is 
denoted  as  follows: 
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[ Hik(t)  =l^lwo«xp  r0j  pxp  lQ]pxw,*  (135) 

This  process  is  called  the  Principal  Component  Analysis.  Theoretically,  this  rank  number  p  is  equal 
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to  twice  of  the  number  of  modes. 


The  minimal  realization  process  begins  from  this  point.  Premultiplying  Eqs.(123)  and  (124)  by  the 
transfer  matrix  [  P  f : 

[/,F[J4(0]«[/,F[*lr«Mi  [M  036) 

[Pr[H*(t+£it)]  =  [Pr  [¥]  feA'A‘i  [L]  (137) 

Notejhat  the  matrix  product  [  P  F  [  ¥  ]  represents  the  reduced  mode  shapes  and  can  be  represented 
by  ( ].  Premultiplying  Eq.(136)  by  the  generalized  inverse  of  [  J 

[Fr  [Pf  [//*(f)]=  re^j  [L]  (138) 

The  right  hand  side  of  this  equation  is  a  common  term  in  Eq.(137).  Substituting  Eq.(138)  into 
Eq.(137)  yields: 

[/,F[/4(f+A/)]  =  [/,Fm  r^'A‘j  K  f[Pf  [#4(0]  (139) 

Simplifying: 

[/4(#+A0]  =  K  ]  r^rA‘j  l*  r  [#4(0]  (140) 

where: 

.  [/4(^Ar)]  =  [Pf  [/4(f+Af)l 

•  [/4(0]  =  [PF[H*(0] 

The  normal  equation  of  Eq.(140)  is: 

[ #4(f  *a t)  1  [ #4(0 F  =  [ ¥ ]  r eArAtj  [¥y[ #4(0 ][ /4(f) f  (i4i) 

To  obtain  the  system  matrix  of  an  eigenvalue  problem,  postmultiply  Eq.(141)  by  the  generalized 

inverse.  [  [#4(0  ]  [#4(0  F  Y 

[/4(f+Af)][/4(oF([f4(0][/4(f)FF  =  [«r ]  feArA‘j  [¥ r  042) 

Postmultiplying  both  sides  of  Eq.  (142)  by  [  ]  and  simplifying  gives : 

[B  )  r eArA*j  =  [tf  ]  re*'j  (143) 

where: 

•  l  B  ]  -  [  /4(f  *a  f )  ]  [  #4(0  f  l  [  #4(0  ]  t  #4(0  F  r 

This  is  an  eigenvalue  problem;  the  size  of  the  system  matrix  [  B  ]  is  now  p  xp  ,  which  is  smaller 
than  N0i  x  N„i  . 

This  is  the  Minimal  Realization  Approach.  The  eigenvectors  are  translated  into  the  mode  shapes  as 
follows: 

1*1  =  (P1K]  (144) 


Transfer  matrix  [  P  ]  is  obtained  from  the  correlation  matrix  by  solving  the  following  eigenvalue 


problem. 


itf*(o][tf*(orm«m  r0*j  045) 

This  algorithm  is  referred  to  as  the  Multiple-Reference  Ibrahim  Time  Domain. 

3.3.3  Modified  Multiple-Reference  Ibrahim  Time  Domain 

In  practice,  the  matrix  size  of  [  HA(t)  ]  [  HA(t)  F  (  N0i  x  Nai  )  is  too  large  to  be  installed  in  many 
mini-computers.  Only  a  few  measurement  data  can  be  included.  This  difficulty  is  overcome  by  some 
matrix  manipulations. 

Transposing  Eq.(123)  and  Eq.(124) 

[**<DF-irF'eMJ  t*r  <146> 

[//*(r+AOF  =  [rr  feAfA‘j  feA,*j  [¥F  (147) 

Decomposing  the  matrix  [  HA(t)  F  as  follows: 

[tf*(01  =  m  fD  j  [Qf  (148) 

or 

[«*(OF-[G]  rOj  [pf  (149) 

Using  this  matrix  [  Q  ]  as  the  transformation  matrix,  Eq.(146)  and  Eq.(147)  can  be  modified  as  fol¬ 
lows: 

[QF[tf*(OF  =  lQF[rF  l*F  (iso) 

or. 

I «i<0 F -  \V7  [*F  (151) 

where: 

.  [//*(D]-l(?Flff*(0F 

•  [L  F  =  [QF(^F  =  reduced  modal  participation  matrix 
Likewise : 

[C?F[J4('+A/)F  =  [(?F[^F  fe*rA‘j  r^j  [»F  (152) 

or, 

[  J&(f  *A  /)  F  =  ( r  F  r eArA< j  r  eMj  [  *  F  (153) 

where: 

.  [ffi(/*AO]-[GF[«*(/*AOF 

Premultiplying  Eq.(151)  by  [  [  L'  F  ]+  : 
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[[l'  rr  1*4(01=  r^Mj  [«r  054) 

The  right  hand  side  of  Eq.(154)  is  a  common  term  in  Eq.(152).  Substituting  Eq.(154)  into  (153)  yields 
[/4(^A0F  =  [I"Fre*'A‘j  [[Z/Fr[/4(0]  055) 

The  normal  equation  of  Eq.(155)  is 

[/4(*+ao ][*&(/) F  =  [i/F  r^rA‘j  [[// Fr[*4(0][/4(0F  056) 

Taking  a  generalized  inverse  of  |[  //*(/)  ]  [  f/»(t)  to  obtain  the  system  matrix  of  an  eigenvalue 
problem: 

[/4(^A0][/4(0F[[/4(0][/4(0Ff  =  [rF  reArA‘j  [[r'rr  (m 

Postmultipiying  Eq.  (157'  by  [  L'  F  and  simplifying 

[B)[U  f  =  [17  F  reA,A‘j 


where: 

•  [  B  ]  =  [  /4(f  +a  t )  ]  [  /4(f)  F  [  [  j&(0  ]  [  /4(f)  Ff 

This  is  the  eigenvalue  problem  to  be  solved;  the  size  of  this  system  matrix  [B  ]  is  pxp  .In  this 
algorithm,  modal  participation  factors  are  calculated  as  eigenvectors.  They  are  transformed  as  fol¬ 
lows: 


[LF  =  [G][jC/F  (158) 

The  transfer  matrix  [  Q  ]  is  calculated  from  the  correlation  matrix  by  the  following  eigenvalue  prob¬ 
lem: 

[«*(0F[^)ne]=[e]rfl5j  d59) 

The  size  of  the  matrix  [  [  //*(/)  F  [  W*(f)  ]J  *s  Af.fcxfyfc  .  This  size  is  usually  smaller  than 

Nai  xN0i  in  Eq.(145),  and  this  algorithm  is,  therefore,  easy  to  implement  on  mini-computers.  This 
algorithm  is  referred  to  as  the  Modified  Multiple-Reference  Ibrahim  Time  Domain.  The  modal  vec¬ 
tors  can  be  calculated  by  using  the  partial  fraction  method  or  by  setting  up  an  overdetermined  system 
using  Eq.  (123).  This  overdetermined  system,  that  has  as  unknowns  only  the  modal  vectors,  can  be 
solved  in  a  least  squares  sense.  These  two  methods  have  been  explained  in  detail  in  the  Polyreference 
Time  Domain  Technique. 

3.3.4  Conclusions 

Two  approaches  of  Multiple-Reference  Ibrahim  Time  Domain  modal  parameter  estimation  tech¬ 
niques,  named  Multiple-Reference  Ibrahim  Time  Domain  and  Modified  Multiple-Reference 
Ibrahim  Time  Domain,  were  developed,  and  attained  sufficient  accuracy  for  modal  parameter  esti¬ 
mation.  As  an  additional  tool  to  predict  the  convergence  of  the  calculated  modal  parameters, 
the  Modal  Confidence  Factor  (MCF)  was  also  included. 


There  is  a  limitation  on  the  developments  of  the  Multiple-Reference  Ibrahim  Time  Domain.  This 
technique  involves  the  formulation  of  a  very  large  system  matrix  for  an  eigenvalue  problem  in 
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order  to  include  alt  measurement  data  Since  several  data  time  shifts  must  be  taken  in  the  system 
matrix  in  order  to  increase  the  rank  of  the  system  matrix,  this  matrix  is  sometimes  larger  than  100 
x  100.  This  limitation  is  avoided  in  the  Modified  Multiple-Reference  Ibrahim  Time  Domain  algo¬ 
rithm.  The  size  of  the  system  matrix  is  always  fixed  to  an  arbitrary  dimension;  usually  this  size  is  lim¬ 
ited  to  40  x  40  which  yields  20  modes. 

The  condition  number  and  rank  chart  is  also  developed  for  the  selection  of  the  number  of  modes  and 
the  effective  rank  in  order  to  better  estimate  the  modal  parameters.  The  rank  chart  usually  indicates 
the  number  of  modes  in  the  measurement  data  within  a  certain  frequency  range,  and  the  rank 
number  at  this  mode  indication  normally  satisfies  the  required  accuracy  for  modal  parameters. 
Where  more  precise  results  are  required,  the  selection  of  a  higher  rank  will  yield  better  results  for 
modal  parameter  estimates  and  convergence  of  MCFs. 

The  Multiple-Reference  Ibrahim  Time  Domain  algorithm  calculates  resonant  frequencies,  damping 
factors  and  mode  shapes  during  the  first  stage  of  modal  parameter  estimation.  The  residues  are 
obtained  with  respect  to  a  specific  reference  point  using  the  normal  equations  method  at  the  second 
stage.  Since  this  process  is  done  in  the  time  domain,  residual  effects  are  not  considered  in  this  algo¬ 
rithm. 

The  Modified  Multiple-Reference  Ibrahim  Time  Domain  obtains  modal  participation  factors  during 
the  first  stage  estimation  the  same  as  the  Polyreference  Time  Domain  technique.  These  modal  parti¬ 
cipation  factors  are  identical  to  those  of  Polyreference  Time  Domain;  therefore,  the  second  stage 
residue  calculation  is  performed  by  a  general  polyreference  residue  calculation  algorithm.  This  algo¬ 
rithm  includes  time  domain  and  frequency  domain  normal  equation  methods.  Section  3.1,  and  also 
has  capability  to  obtain  real  normal  modes.  These  options,  therefore,  are  also  available  with  the 
Modified  Multiple-Reference  Ibrahim  Time  Domain  technique. 

Time  domain  modal  parameter  estimation  algorithms  usually  are  more  numerically  robust  solutions 
than  frequency  domain  algorithms  for  lightly  damped  systems.  The  presented  modal  parameter  esti¬ 
mation  techniques  in  this  section,  Multiple-Reference  Ibrahim  Time  Domain  and  Modified 
Multiple-Reference  Ibrahim  Time  Domain,  both  use  time  domain  data.  For  this  reason,  these 
methods  are  numerically  very  robust  and  consequently  the  Modified  Multiple-Reference  Ibrahim 
Time  Domain  technique  is  comparable  to  the  Polyreference  Time  Domain  technique.  The  Modified 
Multiple-Reference  Ibrahim  Time  Domain  technique  is  able  to  detect  repeated  roots,  since  it  is  a 
multiple-reference  algorithm,  and  has  the  advantage  to  produce  less  computational  poles  in  com¬ 
parison  with  the  Polyreference  Time  Domain  technique. 
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3.4  Multiple-Reference  Orthogonal  Polynomial 


'•"'i 


3.4.1  Introduction 


The  recent  trend  of  experimental  modal  analysis  is  to  apply  a  multiple  excitation  force  and  to  meas¬ 
ure  multiple  response  functions  l51).  This  approach  results  in  a  more  uniform  distribution  of  the  exci¬ 
tation  energy  applied  to  the  structure,  and  hence,  more  uniform  response  level;  thus,  higher  quality 
and  more  consistent  data  can  be  obtained.  Another  practical  advantage  is  that  the  data  acquisition 
time  can  be  reduced  if  the  appropriate  equipment  is  available. 

In  the  last  few  years,  several  "multiple-reference"  parameter  estimation  methods  have  been 
developed  to  take  advantage  of  the  consistency  of  the  frequency  response  measurements  relative  to 
different  reference  locations.  Such  methods  include  Polyreference  Time  and  Frequency  Domain 
Method!18-42),  Eigensystem  Realization  Algorithm!34),  and  Direct  Parameter  Estimation!68).  With 
these  methods,  great  improvement  in  identifying  experimental  modal  models  have  been  achieved. 
Moreover,  the  multiple  reference  approach  can  be  used  to  identify  highly  coupled  and  repeated 
modes,  which  is  difficult  or  even  impossible  to  attack  for  single  reference  approach. 

Most  of  the  early  multiple  reference  work  was  based  upon  time  domain  algorithms.  These  algorithms 
were  characterized  by  numerical  stability,  ease  of  use,  and  were  particularly  good  at  analyzing  wide 
frequency  bands.  However,  there  are  several  important  limitations: 

•  Averaging  is  often  performed  in  the  frequency  domain  and  Fourier  transformed  into  the  time 
domain  for  analysis.  For  this  case,  time  domain  "leakage"  can  be  a  problem,  particularly  for  small 
frequency  bands. 

•  These  time  domain  methods  are  based  upon  sampled  data,  with  the  result  that  all  of  the  com¬ 
puted  frequencies  fall  into  the  Nyquist  band  (from  0  to  1/2  of  the  sampling  frequency).  This  can 
make  it  more  difficult  to  sort  out  computational  modes.  This  is  becoming  increasing  more  impor¬ 
tant  for  the  application  of  active  control  of  large  structures  (Space  Station,  etc)  where 
modal/servo  systems  exist.  These  system  have  both  low  and  high  damped  poles  making  it  particu¬ 
larly  difficult  to  separate  the  computational  poles  from  the  real  poles. 

•  Sine  testing  is  becoming  more  popular  with  the  advent  of  systems  capable  of  making  many  simul¬ 
taneous  measurements.  In  sine  testing,  the  frequency  interval  may  not  be  equal  making  it  difficult 
to  transform  into  the  time  domain. 

•  It  is  difficult  to  weight  the  data  in  a  given  frequency  band  with  time  domain  methods. 

As  a  result,  there  has  been  a  great  deal  of  recent  interest  in  developing  multiple  reference  frequency 
domain  algorithms  as  witnessed  by  the  recent  development  of  the  Frequency  Domain  Polyreference 
and  Eigensystem  Realization  Algorithm  (ERA).  The  rational  fraction  polynomial  is  one  of  the  fre¬ 
quency  domain  algorithms,  which  can  be  extended  to  multiple  references  formulation  with  the  advan¬ 
tage  of  increasing  frequency  resolution  by  the  use  of  multiple  reference  data 

The  general  formulation  for  the  Frequency  Response  Function  (FRF)  is  in  the  form  of  rational  frac¬ 
tion  polynomials.  In  the  frequency  domain  parameter  estimation,  some  numerical  difficulty  prevents 
one  from  using  this  simple  formulation  directly.  One  of  the  difficulties  is  that  the  formulated  equa¬ 
tions  are  usually  ill-conditioned,  resulting  from  the  dynamic  range  of  the  polynomial  terms.  Thus, 
one  of  the  natural  extensions  is  to  use  orthogonal  polynomials  instead  of  ordinary  power  polynomials. 
As  the  property  of  orthogonality  is  exercised,  the  normal  equation  used  in  the  solution  process  is 
easily  decoupled,  and  the  size  of  the  system  matrix  is  greatly  reduced 

In  this  section,  a  rational  fraction  formulation  in  orthogonal  polynomials,  based  on  the  earlier 
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derivation  by  Richardson  and  FormentiP7!,  and  Hout69],  is  extended  to  a  matrix  formulation,  ie.  A 
matrix  Auto-Regressive  Moving  Average  (ARMA)  model  in  the  Laplace  domain!40].  This  method 
can  include  multiple  reference  FRF  calculations  and  thus  is  able  to  detect  highly  coupled  and 
repeated  modes.  The  matrix  polynomial  coefficients  of  the  AR  part  is  first  estimated  in  least  square 
sense,  and  then  the  system  poles  can  be  solved  from  the  characteristic  equation  in  companion  matrix 
form.  The  residue  matrices  can  be  calculated  from  the  matrix  polynomial  coefficients  of  the  MA  part, 
which  is  estimated  in  the  same  least  square  formulation  as  AR  part.  Therefore,  the  global  estimation 
of  both  residues  as  well  as  poles  is  achieved. 

A  singular  value  decomposition  method  has  been  used  to  develop  a  complex  mode  indicator  function 
(CMIF)!70],  which  can  be  used  to  help  determine  the  number  of  poles  before  the  analysis.  The  CMIF 
is  formed  by  performing  a  singular  value  decomposition  of  all  of  the  FRF’s  at  each  spectral  line.  The 
singular  values  are  plotted  as  functions  of  frequency.  These  singular  values  can  also  be  used  to  gen¬ 
erate  appropriate  weighting  functions  in  order  to  enhance  some  certain  frequency  band  and/or  to 
enhance  weak  modes. 


3.4.2  Theory 


For  a  N  degree-of-freedom  linear  system  with  viscous  damping,  one  can  formulate  an  Auto- 
Regressive  Moving  Average  (ARMA)  model  of  order  (m,n)  as 


£[<**]**  [//(*)]=£[&*]** 
k=£)  fc=0 


(160) 


or  in  orthogonal  polynomial  basis 


SKI  I Pk\  [H(s))  =  £  [ft]  \pk] 

*=o  k=a 


(161) 


where: 


•  A  degree-of-freedom  of  the  system  or  the  number  of  modes. 

•  A,-  number  of  excitation  location. 

•  A*  number  of  response  point. 

•  m  order  of  matrix  polynomial  chosen  in  the  Auto-Regressive  part.  mJV(  >  A . 

•  n  order  of  matrix  polynomial  chosen  in  the  Moving- Average  part,  n  >  m  +2 . 

•  [#(*)]  transfer  function  matrix  of  size  A,-  by  A0. 

•  [a*]  matrix  polynomial  coefficient  of  size  by  A,. 

•  [6fc]  matrix  polynomial  coefficient  of  size  A,  by  A„. 

•  [a*]  matrix  polynomial  coefficient  for  orthogonal  polynomials  of  size  A,  by  A,-. 

•  [>9fc]  matrix  polynomial  coefficient  for  orthogonal  polynomials  of  size  A,  by  Na. 

•  \pk]  orthogonal  polynomial  of  order  k,  which  is  a  scalar  function. 


Since  [a*]  and  [f>k]  are  both  unknown  parameters  in  Eq.(161),  one  can  always  choose  one  of  the 
matrix  coefficients  to  be  the  identity  matrix.  By  choosing  [o^,]  =  [/],  Eq.(161)  can  be  arranged  and 
rewritten  in  matrix  form  for  all  frequency  channels. 
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(162) 


where: 


•  [R*F  product  of  pk  and  identity  matrix  [/]  of  size  by  N{. 

•  PVF  product  of  pk  and  identity  matrix  [/]  of  size  N0by  N0. 


In  order  to  solve  for  the  matrix  coefficients,  one  can  formulate  the  normal  equation  and  further 
decouple  it  by  taking  the  advantage  of  the  orthogonal  property. 


!3  Iff]  [HI- 


mi 


(163) 


or 


([Xf\X]-[Y])[a]  =  [Xf  K,]-[yj 

M  -KJ-OTH 


(164) 

(165) 


where: 


•  [X]  =  Re 


-  EFoWFol*  ...  - 


•=i 
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*  Y,[P*W\B[Po\H 

i=i 


•  [*,]  =  Re 


E  P*o  ‘f[Hf[Pmf 

•=i 


i=i 


•  [K]  =  Re 


e raw] [h]b[p0]s  ...  z[Pom[HriPm*r 

*=i  ■=! 


e  [nnp0]H...  e  [n]B[pm.i\B 
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•  [Yn]  =  Re 


[Oof 

Vof 

•  N  = 

K- if 

L9]  = 

[finf 

Note  here  that  the  size  of  the  system  matrix  in  Eq.(164)  is  mN{  by  mN{ ,  which  is  significantly  smaller 
than  the  size  of  the  previous  system  matrix  in  Eq.(163). 

The  orthogonal  polynomial  with  any  weighting  function  q(s)  is  defined  and  normalized  as: 

E P&iUbdPktof  -  6*  (166) 

i=l 

A  step  by  step  procedure  for  generating  the  orthogonal  polynomials  is  developed  by  Forsythet73).  A 
simplified  Forsythe  orthogonal  polynomials,  namely  Real  Half-Function  Orthogonal  Polynomials)37), 
for  complex  variable  function  can  be  derived  as: 


Normalized  function 
5-i(h)  =  0 


S0(>W)  =  -r 

«o 


Nonnormalized  function 
R.  i(h)  =  0 
Ro(h)  =  1 


5x(h) 


53(h)  = 


<*o  =  (2  E*)s 
i=i 

Ri(h) 


^i=(2E^iU)aft)3 

1=1 

r2(h) 


Ri(h)  =  H-So(“v) 

V0  -  0 

R3(h)=  h  51(h)*»/i  5o(h) 
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A  transformation  matrix  [r]  can  also  be  calculated  in  the  generating  procedure.  Thus  the  power 
polynomial  coefficients  can  easily  be  calculated  from  orthogonal  polynomial  coefficients  [a]  and 
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where,  the  subscript  ik  indicate  only  one  of  the  entiy  ik  of  the  matrix  polynomial  coefficients  is 
considered. 
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Once  the  [a*]’s  are  found,  one  can  solve  for  the  poles,  natural  frequencies  and  dampings  characteris¬ 
tics,  by  setting  the  determinant  of  the  characteristic  equation  equal  to  zero. 


E  =  0  (169) 

frd> 

Instead  of  solving  the  matrix  characteristic  equation,  one  can  reformulate  Eq.(169)  in  the  companion 
matrix  form  and  solve  the  eigenvalue  problem  for  the  poles. 
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Note  that  in  this  formulation  {v}  is  one  of  the  columns  of  the  modal  participation  matrix,  which  is 
proportional  to  the  modal  coefficients  at  the  reference  degrees  of  freedom  for  one  mode. 

Assuming  that  no  exact  repeated  roots  exist,  which  is  always  the  case,  then  the  residue  matrix  for 
mode  r  can  be  defined  as 

lim  (s-Ar)[//(s)]  (171) 

by  applying  Eq.(160) 

t/U*  lim  (s  -  A  r)  [a  (s)J*1  [b(s))  (172) 

•-•A, 

where: 

•  [a(i)]=  £  fa*!** 


ksaO 


Since  [a(s)]  is  singular  as  s-»Ar  ,  one  can  formulate  a  Singular  Value  Decomposition  to  determine 
the  generalize  inverse  of  [a  (j  )] 

[a  (s)] -[(/(*)]  [E^)]^)]*  (173) 

or 

[«  (•*)]  *  E  *.(*)  {«<(■»)}  {  v.(s)}*  (174) 

isl 


where: 


ms)] 


*l(s) 


<r3(s) 


°nSs) 


singular  value  matrix,  which  is  a  diagonal  matrix. 


[U(s)]  -  J{Ui(s)}  {ua(s)}  •  •  •  {wW|(s)}J  left  singular  matrix,  which  is  a  unitary  matrix. 
[//(*)]*  [U(s)]  =  [/] 

[^(s)]  =  [{Vi(s)}  (v2(s)}  •  •  •  {vW((s)}j  right  singular  matrix,  which  is  a  unitary  matrix. 


-63- 


•[K(5)]S[K(5)]  =  [/] 


nvwTiiy  jww  *r  w.nn ji#\ jv-h  ^  ’  «'<m>vir 


The  generalized  inverse  is  obtained  as 

[am1  -[vwmsn1  [u(s))a 


(175) 


or 


w. 
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H 


(176) 


Applying  Eq.  (176)  to  the  indeterminate  term  in  Eq.(172): 


M 


lim  (s -A r)[a (s)]*1  =  lim  (s-Ar)  ^.(s)1  {v.(s)}{Ki(s)}ir 

«-«Ar  «-**r  ,_1 


(177) 


If  Eq.(174)  is  evaluated  at  s  =  A  r ,  at  least  one  of  the  singular  values,  cr,(A  r) ,  must  be  zero,  which 
indicates  that  the  matrix  ( a  (A ,)  ]  is  not  of  full  rank.  This  is  consistent  with  the  inability  to  deter¬ 
mine  the  eigenvector  in  an  absolute  sense.  For  the  case  of  repeated  roots  at  s  =  A , ,  more  than  one 
of  the  singular  values  will  be  zero.  For  the  case  of  no  repeated  roots,  only  one  of  the  singular  values, 
<r,(Af) ,  will  be  zero.  Thus,  Eq.(177)  reduces  to: 


fs-A  ) 

lim  (j  -A  r)[a  (i)]4  =  lim  ~  K(A,)}{«,(Ar)}g 

#-*Ar  <Ti[S ) 


(178) 


For  the  identity  that  ^(A,)  =  0  ,  the  reciprocal  of  the  indeterminate  term  in  Eq.(178)  can  be 
expressed  as: 


o.(s)  ..  ff,(s)  -ff,(Ar)  d  <7,(s) 

lim  — r~  -  hm 


1*4,  s  -A  P  •— *Ar  J  -  A  r 

Therefore,  Eq.(178)  becomes: 

*_,<**(*) 


ds 
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(179) 


lim  (s-Af)[n(j)r  =  (-  , 

•—a,  as 


-OM^A,)}^,)}* 


(180) 
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In  order  to  determine  — - —  |  ,_*r ,  one  can  take  the  derivative  of  Eq.(174)  as 


d\a  (*)] 


ds 
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ds 

*  ,  d*i(s) 


i=l 


ds 


,  ,  vw  ,  x  </{u.(5)}  ,  ^Vljr  #w  ,  „  ^{v,(s)}fl  x  , 

{UiisMViis)}*  *  ffi(s ) — - - (v,(s)}g  +  <r<(s  ){«,($)} - X - >  I  *=v 


ds 


(181) 


Then,  by  pre-multiplying  Eq.(181)  by  {w,(Ar)}g  ,  post-multiplying  by  {v,(Ar)} ,  and  using  the  ortho¬ 
gonality  property  of  a  unitary  matrix,  only  one  term  remains  on  the  right  hand  side. 


d<Ti(s) 


{«.( Ar)}'^M|..xr{v,(Ar)}-  ^ 


I  »=Xr 


(182) 


Substituting  this  result  into  Eq.(180), 
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(v.(A  ,)}{«.(>,»* 

lun  (j-AJfafs)]-1  - - —■  - - 

{«.(A,)}*  |#sAr  (v,(Ar)} 


Therefore,  the  residue  matrix  in  Eq.(172)  can  be  calculated  as, 

KJ“ - .r  T-Ti -  [0(Af)J 

{u.(Xr)}B  dl2££\  {v.(Xf)} 


3.4.3  Conclusions 


The  Multiple-Reference  Orthogonal  Polynomial  has  demonstrated  some  attractive  characteristics  for 
modal  parameter  estimation.  Some  of  these  are: 

•  It  can  use  any  frequency  spacing  data,  and  thus  is  suitable  for  step  sine  measurement  data 

•  Works  well  in  narrow  frequency  bands 

•  Global  parameter  estimation  for  all  FRF  measurement. 

•  It  includes  multiple  reference  measurements  calculation  and,  therefore,  it  is  good  for  system  with 
highly  coupled  and  repeated  modes 

•  Easy  to  increase  the  order  of  the  polynomials  for  residual  calculation  and  error  compensation. 

•  Weighting  as  a  function  of  frequency  can  be  used 

•  Additional  weighting  function  can  be  used  to  detect  or  enhance  weak  modes. 

•  By  direct  use  of  FRF  data,  time  domain  "leakage"  is  avoided. 

•  Fewer  computational  modes  introduced,  thus  it  is  suitable  for  combined  modal/servo  systems 
which  can  include  both  low  and  high  damping  modes. 

•  CMIF  can  be  used  to  indicate  close  or  repeated  eigenvalues  before  the  parameter  estimation  pro¬ 
cedure,  so  that  the  proper  order  of  the  polynomials  can  be  chosen. 


3.5  Multi-Mac 


3.5.1  Introduction 

Multi-Mac  is  a  spatial  domain  method  of  determining  modal  parameters  based  on  multiple-reference 
frequency  response  functions.  Multi-Mac  is  an  extension  of  the  concept  of  Modal  Assurance  Cri¬ 
terion  (MAC).  MAC  is  a  calculation  that  is  used  to  gain  confidence  in  the  estimates  of  modal  vec¬ 
tors  either  by  verifying  that  the  estimates  of  different  modal  vectors  are  unique,  that  normal  modes 
have  been  estimated,  or  that  estimates  from  different  rows  or  columns  of  the  residue  matrix  for  the 
same  mode  are  identical.  Multi-Mac  uses  principal  component  analysis  to  determine  the  number  of 
independent  vectors  that  make  up  the  residue  vectors  calculated  from  a  number  of  reference  loca¬ 
tions.  The  number  of  significant  eigenvalues  is  used  to  determine  the  number  of  independent  vectors. 
Using  the  eigenvectors  from  the  principal  component  analysis,  a  set  of  orthogonal  mode  shapes  can 
be  calculated. 

3.5.2  Theory 

For  linear  systems  with  normal  modes,  the  residue  of  any  input/output  combination  is  made  up  of 
three  parts,  the  modal  vector  at  the  input  location,  the  modal  vector  at  the  output  location,  and  the 
scaling  constant  for  that  model50-54-76"78). 

Siptr  =  Qr\VW  (185) 

Where: 

•  Qr  =  Scaling  constant  for  mode  r 

•  i>f,  s  Modal  vector  for  location  p  of  mode  r 

•  W  =  Modal  vector  for  location  q  of  mode  r 

•  A,,,, *»  Residue 

Equation  (185)  can  be  expanded  to  form  the  complete  residue  matrix  for  any  mode  r. 

'kr'lhr  - 

iM*  tMbr  - 

[Ar)  =  Qr  086) 

VWtkr  VwVbr  VwVbr  - 

From  Equations  (186)  it  can  be  seen  that  there  must  be  a  linear  relationship  between  any  row  or  any 
column  of  the  residue  matrix  for  a  particular  mode.  For  example,  if  the  structure  was  excited  at 
points  1  and  2,  then  in  column  one  every  residue  would  have  in  common.  Likewise,  every  residue 
in  row  two  would  have  in  common.  The  modal  vectors  at  the  output  loactions,  through  Vw- 
must  therefore  be  related  since  they  define  the  same  mode.  If  this  linear  relationship  does  not  exist, 
either  measurement  errors  have  contaminated  the  data,  the  modes  are  closely  coupled  or,  there  is  a 
repeated  root  at  that  pole  location  so  that  the  residue  is  a  linear  combination  based  on  input  loca¬ 
tion. 

Also,  from  the  basic  theory  of  modal  analysis,  the  modal  vectors  of  two  different  modes  must  be 
orthogonal  with  each  other  and  a  weighting  matrix  such  as  the  mass  matrix. 


If  estimates  of  the  residues  for  the  same  mode  but  for  different  rows  or  columns  is  used,  a  principal 
component  analysis  can  be  used  to  determine  the  number  of  independent  vectors  that  make  up  those 
residues.  If  one  significant  eigenvalue  is  found  from  the  principal  component  analysis,  then  one 
mode  is  presented.  If  more  than  one  significant  eigenvalue  is  calculated,  then  there  are  more  than 
one  set  of  independent  vectors  that  make  up  the  residue  matrix  for  that  mode. 

Equation  (186)  can  be  used  to  set  up  the  principal  component  analysis.  The  size  of  the  matrices  is 
determined  by  the  number  of  modal  vector  estimates  and  the  number  of  measurement  locations. 
The  residue  matrix  can  have  more  than  one  estimate  of  the  modal  vector  for  the  same  input/output 
locations.  These  estimates  could  be  from  different  curve  fitting  algorithms.  The  matrix  [U]  is  the 
matrix  of  eigenvectors  for  each  eigenvalue. 

\4r]  \Ar]a  =  [W]  [U]  \  Aj  [Uf  [Wf  (187) 

Where: 

•  [A]  =  Residue  matrix 

•  [W]  =  Diagonal  weighting  matrix 

•  [U]  =  Eigenvector  matrix 

•  \  AJ  =  Diagonal  eigenvalue  matrix 

Equation  (187)  will  yield  as  many  eigenvalues  as  the  number  of  rows  in  the  residue  matrix.  Some  of 
these  eigenvalues  will  be  zero  or  of  insignificant  value.  The  number  of  significant  eigenvalues  will 
indicate  the  number  of  independent  vectors  that  make  up  the  residues  at  that  frequency.  Associated 
with  each  eigenvalue  will  be  the  eigenvectors  for  that  eigenvalue  which  will  number  as  many  as  the 
number  of  rows  in  the  residue  matrix.  The  weighting  matrix  that  was  used  in  Eq.  (187)  for  this  work 
was  the  identity  matrix.  A  better  weighting  matrix  would  be  to  use  the  inverse  of  the  square  root  of 
the  mass  matrix  or  an  estimate  of  the  inverse  of  the  square  root  of  the  mass  matrix  from  a  finite  ele¬ 
ment  analysis.  The  matrix  could  also  be  an  error  matrix  or  a  matrix  that  would  allow  different  types 
of  data  in  the  residue  matrix.  For  example,  if  both  acceleration  and  displacement  data  was  contained 
in  the  residue  matrix,  the  multiplication  of  acceleration  and  displacement  would  not  be  dimensionally 
correct.  But,  the  weighting  matrix  could  be  used  to  make  the  multiplication  dimensionally  correct!70]. 

Any  estimate  of  the  residue  may  be  used  in  this  procedure.  Since,  at  resonance,  the  quadrature  part 
of  the  frequency  response  function  is  proportional  to  the  residue,  the  simplest  method  is  to  use  the 
peak  of  the  imaginary  part  of  the  frequency  response  estimate  as  an  estimate  of  the  residue!16!.  If 
leakage  is  present  in  the  data,  it  is  advantageous  to  use  spectral  lines  that  are  adjacent  to  the  highest 
peak  in  the  imaginary  part.  These  spectral  lines  will  be  less  contaminated  by  leakage.  Therefore,  the 
multi-mac  calculation  can  be  used  before  any  sophisticated  parameter  estimation  is  attempted  to 
define  the  number  of  independent  vectors  at  each  pole  location  or  can  be  used  after  all  parameter 
estimation  has  been  completed  to  gain  confidence  in  the  estimate  modal  vectors. 

Once  the  eigenvalues  and  eigenvectors  have  been  calculated  by  the  principal  component  analysis,  the 
residues  can  be  transformed  to  a  new  coordinate  space  which  are  mathematically  guaranteed  to  be 
orthogonal  with  each  other  and  the  assumed  weighting  matrix.  Therefore,  if  more  than  one  eigen¬ 
value  of  significant  value  is  found,  the  orthogonal  mode  shapes  that  make  up  the  estimate  of  the  resi¬ 
due  at  that  frequency  can  be  computed.  For  example,  if  there  are  repeated  roots,  the  residues  that 
make  up  that  repeated  root  will  be  calculated.  If  the  modes  are  heavily  coupled,  unique  modes  will 
be  estimated.  But,  the  transformed  modes  will  have  the  contamination  of  other  modes  removed  and 
will  be  othogonal  with  each  other. 

Equation  (187)  defines  the  transformation  from  the  original  residues  to  the  new  space.  This 
transformation  would  be  used  for  each  eigenvalue  of  significant  value.  Therefore,  if  more  than  one 


eigenvalue  of  significant  value  is  found,  the  transformation  will  yield  orthogonal  estimates  of  the  resi¬ 
dues.  If  more  than  one  estimate  of  the  residue  is  used  for  each  row  or  column,  these  transformed 
estimates  are  summed  together  to  reduce  the  variance.  The  matrix  multiplication  will  therefore  yield 
one  residue  for  each  measurement  location.  The  variance  of  the  modal  vectors  will  be  reduced  in  a 
least  squares  sense  by  [1  In  ]^3. 


(188) 

3.5.3  Experimental  Procedure 

In  practice,  a  number  of  reference  locations  are  used  to  estimate  the  frequency  response  functions. 
These  can  be  from  either  a  multiple  input  test  or  an  impact  test  with  a  number  of  reference 
accelerometers.  In  this  way,  a  number  of  rows  equal  to  the  number  of  reference  accelerometers  or  a 
number  of  columns  equal  to  the  number  of  inputs  of  the  residue  matrix  will  be  estimated.  To  get 
good  estimates  of  the  natural  frequencies  of  the  system,  a  summation  of  the  power  spectrum  of  the 
quadrature  response  is  calculated!16].  This  summation,  which  can  use  all  frequency  response  estimates 
or  a  subset  of  the  frequency  response  estimates,  enhances  the  global  modes  of  the  system  and  give 
less  weight  to  any  local  modes. 

From  these  starting  values,  the  quadrature  part  of  the  frequency  response  estimates  are  used.  In 
practice,  the  peak  location  and  one  or  two  spectral  lines  on  either  side  of  the  peak  is  used.  In  this 
way,  the  errors  caused  by  leakage  can  be  reduced.  If  there  is  no  leakage  or  if  a  large  number  of  refer¬ 
ences  is  used,  it  is  not  necessary  to  use  more  than  the  peak  value.  These  three  or  five  frequencies  are 
used  as  an  estimate  of  the  residue  for  that  mode  for  a  given  reference  location. 

A  principal  component  analysis,  which  is  the  number  of  spectral  lines  times  the  number  of  references 
in  size,  is  then  done  on  these  estimates  of  the  residue.  This  will  yield  a  number  of  eigenvalues  equal 
to  the  number  of  spectral  lines  times  the  number  of  references.  If  only  one  eigenvalue  of  significant 
value  is  calculated,  then  there  is  one  independent  vector  that  makes  up  the  residue  vectors  at  those 
frequencies.  If  more  that  one  significant  eigenvalue  is  found  then  there  are  that  many  independent 
vectors  that  make  up  the  residue  vectors  at  that  frequency. 

The  residue  vectors  can  then  be  transformed  by  the  eigenvectors  found  from  the  principal  component 
analysis.  These  new  mode  shapes  are  mathematically  guaranteed  to  be  orthogonal  with  each  other. 
These  transformed  residues  could  be  weighted  by  the  mass  matrix  to  yield  mode  shapes  of  the  system. 

Using  these  transformed  estimates  of  the  residue  vectors,  an  enhanced  frequency  response  function 
can  be  computed  which  can  be  fit  for  estimates  of  the  frequency  and  damping  of  that  model50!. 
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4.  CASE  STUDIES 


4.1  Circular  Plate 

Utilizing  impact  testing,  data  was  taken  on  a  circular  plate  structure.  Seven  accelerometers  were  per¬ 
manently  mounted  on  the  structure,  while  the  force  was  moved  to  36  different  locations.  Therefore, 
a  data  set  was  obtained  with  seven  references.  Figure  5  shows  the  circular  plate  modal  model  with 
the  seven  references  points  indicated. 
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Figure  5.  Circular  Plate  Model 

The  following  algorithms  were  used  to  analyze  the  data  for  frequency  and  damping  characteristics: 

•  Polyreference  Time  Domain  (PTD) 

•  Polyreference  Frequency  Domain  (PFD) 

•  Orthogonal  Polynomial  (OP) 

•  Multiple-Reference  Ibrahim  Time  Domain  (MRITD) 

•  Modified  Multiple-Reference  Ibrahim  Time  Domain  (MMRITD) 

Due  to  memory  limitations  in  the  implementation  of  the  algorithms,  the  data  was  analyzed  with  only 
six  references.  The  references  used  for  the  analysis  were  arbitrarily  chosen.  The  data  was  obtained 
using  baseband  analysis  from  0  Hz.  to  2560  Hz.,  with  a  time  domain  blocksize  of  1024.  Therefore,  the 
frequency  resolution.  A/,  was  5  Hz.  A  typical  driving  point  frequency  response  function  is  shown  in 
Figure  6. 


The  data  was  processed  in  several  analyses.  First,  the  data  was  analyzed  from  spectral  line  50  to  spec¬ 
tral  line  178  or,  from  250  Hz.  to  890  Hz.  The  second  analysis  band  was  from  spectral  line  200  to  spec- 
tal  line  328  or,  from  1000  Hz.  to  1640  Hz.  The  frequency  and  damping  results  of  these  analyses  are 
presented  in  Table  1. 
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TABLE  1.  Circular  Plate  Frequency  and  Damping  Analysis 


Mode 

Frequency/Damping  Estimation  Method 

PTD 

PFD 

OP 

MRITD 

MMRITD 

1 

Freq.  (Hz) 
Damp.  (%) 

362.276 

0.8665 

362.368 

0.8296 

362.474 

0.8876 

362.208 

0.8580 

362.254 

0.868 

2 

Freq.  (Hz) 
Damp.  (%) 

363.612 

0.9368 

363.562 

0.8850 

363.603 

0.9554 

363.629 

0.9013 

363.567 

0.946 

3 

Freq.  (Hz) 
Damp.  (%) 

557.041 

0.5129 

557.124 

0.5683 

557.078 

0.5389 

557.085 

03154 

557.045 

0.523 

4 

761.143 

0.6654 

761.309 

0.6287 

761.185 

0.6667 

761.228 

0.6653 

761.171 

0.663 

5 

Freq.  (Hz) 
Damp.  (%) 

764.138 

0.3371 

764.166 

0.3259 

764.146 

0.3338 

764340 

03335 

764.190 

0.339 

6 

Freq.  (Hz) 
Damp.  (%) 

1223.021 

0.3337 

1223.119 

0.3445 

1223.006 

03357 

1223.033 

03372 

1223.006 

0.334 

D 

Freq.  (Hz) 
Damp.  (%) 

1224.126 

0.3262 

1223.814 

0.3195 

1224.145 

03211 

1224.034 

0.3314 

1224.092 

0332 

8 

Freq.  (Hz) 
Damp.  (%) 

1328.063 

0.5031 

1328.213 

0.5207 

1328.030 

0.4999 

1328.070 

0.4936 

1328.069 

0.496 

9 

Freq.  (Hz) 
Damp.  (%) 

1328.760 

0.4129 

1328.687 

0.4079 

1328.754 

0.4084 

1328.772 

0.4032 

1328.781 

0.408 

Notes: 


PTD  Polyreference  Time  Domain 

PFD  Polyreference  Frequency  Domain 

OP  Orthogonal  Polynomial 

MRlTD  Multiple-Reference  Ibrahim  Time  Domain 

MMRITD  Modified  Multiple-Reference  Ibrahim  Time  Domain 


42  Body-In-White 

Frequency  response  functions  were  obtained  from  a  "body-in-white"  automobile  structure  using 
dual-shaker  random  excitation  with  triaxial  responses.  The  data  set  contains  432  frequency  response 
functions;  from  the  72  measurement  locations,  measured  in  the  three  global  directions,  and  the  two 
references.  Figure  7  shows  the  "body-in-white"  modal  model  with  the  two  reference  points  indicated. 
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Figure  7.  Body-In-' White  Model 

The  frequency  response  functions  were  obtained  using  baseband  analysis  from  0  Hz.  to  60  Hz.,  with  a 
time  domain  blocksize  of  1024.  Therefore,  the  frequency  resolution.  A/,  was  0.117  Hz.  A  typical 
driving  point  frequency  response  function  is  shown  in  Fig.  8. 


-72- 


V 


BgoBBBiimnicmnp^BPiioBiTOBgtTBt  m CTTOTProt  rorenonorora?  ra  pj;  to  ».v 


Figure  8.  Body-In- White  Driving  Point  Frequency  Response  Function 


The  following  algorithms  were  used  to  analyze  the  data  for  frequency  and  damping  characteristics: 

•  Polyreference  Time  Domain  (PTD) 

•  Eigensystem  Realization  Algorithm  (ERA) 

•  Modified  Multiple-Reference  Ibrahim  Time  Domain  (MMRITD) 

The  frequency  and  damping  evaluation  presented  in  Table  2  was  performed  as  a  class  project  at  the 
University  of  Cincinnati!65-79]. 

4.3  Multi-Mac  Experimental  Examples 

To  demonstrate  the  value  of  the  Multi-Mac  calculation,  a  number  of  test  cases  were  evaluated. 

First,  a  H-frame  was  tested.  Figure  9  shows  the  H-frame.  Twelve  reference  accelerometers  were 
used,  four  tri-axial  accelerometers  at  the  four  ends  as  shown  in  Figure  9.  The  H-frame  was  impacted 
at  29  points  in  the  vertical  direction  and  26  points  in  the  horizontal  direction. 

Figure  10  shows  a  typical  driving  point  frequency  response  function  estimate  in  the  range  of  430 
Hertz  to  630  Hertz.  There  appears  to  be  3  or  4  natural  frequencies  in  this  range.  Figure  11  is  a  sum¬ 
mation  of  the  power  spectrum  of  the  quadrature  responses  using  points  1Z,  13Z,  17Z,  29Z,  1Y,  and 
13Y  as  the  references.  This  function  shows  that,  using  these  references,  there  are  at  least  7  global 
modes  in  this  range. 
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TABLE  2.  Body-In* White  Frequency  and  Damping  Analysis 


Mode 

Frequency/Damping  Estimation  Method 

PTD 

ERA 

MMRITD 

1 

Freq.  (Hz) 

23.704 

23.740 

23.693 

Damp.  (%) 

0.6848 

0.9602 

0.637 

Freq.  (Hz) 

32.174 

32.173 

32.175 

z 

Damp.  (%) 

0.6265 

0.6161 

0.625 

Freq.  (Hz) 

35.293 

35.294 

35.293 

3 

Damp.  (%) 

0.4460 

0.4653 

0.445 

Freq.  (Hz) 

40.636 

40.658 

40.663 

4 

Damp.  (%) 

0.5425 

0.4779 

0301 

Freq.  (Hz) 

41.691 

41.713 

41.719 

5 

Damp.  (%) 

0.4263 

03306 

0.404 

Freq.  (Hz) 

42.760 

42.767 

42.750 

6 

Damp.  (%) 

0.3728 

03762 

0374 

Freq.  (Hz) 

43383 

43.382 

43.385 

7 

Damp.  (%) 

0.4357 

0.4264 

0.439 

8 

Freq.  (Hz) 

46.748 

46.784 

46.728 

Damp.  (%) 

0.7213 

0.9369 

0.656 

Freq.  (Hz) 

48.819 

48.955 

48.943 

y 

Damp.  (%) 

0.8949 

0.8937 

0.955 

10 

Freq.  (Hz) 
Damp.  (%) 

53.926 

03771 

53.938 

03487 

n/a 

Notes:  PTD 

Polyreference  Time  Domain 

ERA 

Eigensystem  Realization  Algorithm 

MMRITD 

Modified  Multiple-Reference  Ibrahim  Time  Domain 

n/a 

not  analyzed 

To  demonstrate  the  value  of  Multi-Mac,  a  principal  component  analysis  was  performed  on  each 
mode.  Tables  3  through  9  are  the  first  20  eigenvalues  from  the  principal  component  analysis.  These 
eigenvalues  are  from  five  spectral  lines  and  six  references  locations  using  the  quadrature  part  of  the 
frequency  response  estimate.  The  tables  show  the  most  significant  eigenvalues  as  a  function  of  the 
rate  of  change  in  the  eigenvalues. 

Figures  12  through  25  are  the  transformed  mode  shapes  using  the  eigenvectors  found  from  the  prin¬ 
cipal  component  analysis.  If  more  than  one  significant  eigenvalue  was  found,  then  that  many  modes 
were  calculated.  These  modes  are  orthogonal  with  each  other  due  to  the  transformation.  Although 
none  of  the  modes  are  repeated  roots,  the  Multi-Mac  calculation  shows  the  high  degree  of  modal 
coupling  for  many  of  the  modes  in  this  range.  As  an  example,  the  estimate  of  the  residue  vector  at 
576  Hertz  is  contaminated  by  both  the  mode  at  567  Hertz  and  the  mode  at  581  Hertz.  This  can  be 
seen  by  noting  that  the  one  mode  computed  for  576  Hertz  is  the  same  as  the  567  Hertz  and  the 
modes  at  581  are  also  similar  to  the  576  Hertz  estimate.  But,  as  shown  by  Table  10,  the  estimates  are 
orthogonal  to  each  other.  This  is  true  for  any  of  the  frequencies  that  have  more  than  one  significant 
eigenvalue. 
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TABLE  3.  H-Crame  Eigenvalues  for  479  Hertz 

********************************** 

***  CRITERIA  TO  JUDGE  RANK  *** 
********************************** 


N0.= 

1 

E.VAL- 

• 127152E+03 @*************************** 

N0.= 

2 

E.VAL- 

. 257414E+00@* 

NO.** 

3 

E.VAL- 

. 823126E“01@* 

NO.- 

4 

E.VAL- 

. 314892E-02@* 

NO.- 

5 

E.VAL- 

. 9O9428E-O20* 

NO.» 

6 

E.VAL- 

.39144OE-O20* 

NO.- 

7 

E.VAL- 

. 544588E-O30* 

NO.- 

8 

E.VAL- 

.486O9OE-O30* 

NO .  — 

9 

E.VAL- 

. 318O55E-O30* 

NO .  — 

10 

E. VAL— 

. 138132E-036* 

NO.- 

11 

E.VAL- 

. 930068E-04@* 

NO.  — 

12 

E.VAL— 

. 522677E-O40* 

NO.  — 

13 

E.VAL- 

. 409133E-046* 

NO .  — 

14 

E.VAL— 

. 146818E-04@* 

NO.- 

15 

E.VAL- 

. 102797E-04®* 

NO.* 

16 

E.VAL- 

. 869041E-05@* 

NO.- 

17 

E.VAL— 

.  786539E-05®* 

NO.— 

18 

E.VAL- 

. 689718E-05@* 

NO.— 

19 

E. VAL— 

.  483071E-05@* 

NO.— 

20 

E.VAL- 

.  347911E-05@* 

SELECTED 

RANK  OF 

DD  MATRIX-  1 

a 


Figure  12.  H-frame  Mode  Shape  at  479  Hertz 
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TABLE  4.  H-frame  Eigenvalues  for  492  Hertz 


********************************** 

***  CRITERIA  TO  JUDGE  RANK  *** 
********************************** 


NO.* 

1 

E.VAL* 

. 220459E+02@*** ************************ 

NO.* 

2 

E.VAL* 

. 647909E+00@ ***************** 

NO.* 

3 

E.VAL* 

. 282519E-O10** 

NO.* 

4 

E.VAL* 

. 9619OOE-O20* 

NO.* 

5 

E.VAL* 

.528314E-026*** 

NO.* 

6 

E.VAL* 

. 126037E-02@* 

NO.* 

7 

E.VAL* 

. 843O33E-O30* 

NO.* 

8 

E.VAL* 

. 344217E-036* 

NO.* 

9 

E.VAL* 

. 263810E-03@* 

NO.* 

10 

E.VAL* 

. 162976E-O30* 

NO.* 

11 

E.VAL* 

. 920806E-04@* 

NO.* 

12 

E.VAL* 

. 729814E-04@* 

NO.* 

13 

E.VAL* 

. 427178E-O40* 

NO.* 

14 

E.VAL* 

. 314986E-04@* 

NO.* 

15 

E.VAL* 

. 2O8977E-O40* 

NO.* 

16 

E.VAL* 

. 831395E-05@* 

NO.* 

17 

E.VAL* 

. 512812E-056* 

NO.* 

18 

E.VAL* 

. 376602E-05§* 

NO.* 

19 

E.VAL* 

. 180069E-05@* 

NO.* 

20 

E.VAL* 

. 140457E-05§* 

SELECTED 

RANK  OF 

DD  MATRIX-  2 

Figure  13.  H-frame  First  Mode  Shape  at  492  Hertz 


TABLE  5.  H-Crame  Eigenvalues  for  500  Hertz 


********************************** 

***  CRITERIA  TO  JUDGE  RANK  *** 
********************************** 

. 194 8 3 3E+02@ *********** ****************§ 
. 211252E+OO0*  @ 

. 995534E-O10*  § 

. 158864E-O10*  0 

. 743O9OE-O20*** 

. 518315E-O30* 

. 4OO569E-O30* 

. 341622E-O30* 

. 156928E-O30* 

. 116135E-O30* 

. 68537 6E-O40* 

. 428652E-O40* 

. 2562O7E-O40* 

. 1OO953E-O40* 

. 558249E-O50* 

.469O59E-O50* 

. 313182E-O50* 

. 173O96E-O50* 

. 123O87E-O50* 

. 883329E-O60* 

MATRIX*  1 


NO.* 

1 

E.VAL* 

NO.* 

2 

E.VAL* 

NO.* 

3 

E.VAL* 

NO.* 

4 

E.VAL* 

NO.* 

5 

E.VAL* 

NO.* 

6 

E.VAL* 

NO.* 

7 

E.VAL* 

NO.* 

8 

E.VAL* 

NO.* 

9 

E.VAL* 

NO.* 

10 

E.VAL* 

NO.* 

11 

E.VAL* 

NO.* 

12 

E.VAL* 

NO.* 

13 

E.VAL* 

NO.* 

14 

E.VAL* 

NO.* 

15 

E.VAL* 

NO.* 

16 

E.VAL* 

NO.* 

17 

E.VAL* 

NO.* 

18 

E.VAL* 

NO.* 

19 

E.VAL* 

NO.* 

20 

E.VAL* 

SELECTED 

RANK  OF 

DD 
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Figure  15.  H-frame  Mode  Shape  at  500  Hertz 


TABLE  6.  H-frame  Eigenvalues  for  567  Hertz 


********************************** 

***  CRITERIA  TO  JUDGE  RANK  *** 
********************************** 


NO.* 

1 

E.VAL* 

. 83562 1E+O20*************************** 

NO.* 

2 

E.VAL* 

. 712774E+00@* 

NO.* 

3 

E.VAL* 

. 295527E+006* 

NO.* 

4 

E.VAL* 

. 151509E+00@* 

NO.* 

5 

E.VAL* 

. 336923E-01@* 

NO.* 

6 

E.VAL* 

. 132355E-01§* 

NO.* 

7 

E.VAL* 

. 332226E-02@* 

NO.* 

8 

E.VAL* 

. 2375O6E-O20* 

NO.* 

9 

E.VAL* 

. 11O479E-O20* 

NO.* 

10 

E.VAL* 

. 488129E-O30* 

NO.* 

11 

E.VAL* 

. 400007E-03@* 

NO.* 

12 

E.VAL* 

. 157472E-03@* 

NO.* 

13 

E.VAL* 

. 994866E-O40* 

NO.* 

14 

E.VAL* 

. 5 5080 IE-040* 

NO.* 

15 

E.VAL* 

. 318778E-O40* 

NO.* 

16 

E.VAL* 

. 11143  3E-04  @* 

NO.* 

17 

E.VAL* 

.965316E-O50* 

NO.* 

18 

E.VAL* 

. 7598O9E-O50* 

NO.* 

19 

E.VAL* 

. 52O88QE-O50* 

NO.* 

20 

E.VAL* 

. 3949O4E-O50* 
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Figure  16.  H-frame  Mode  Shape  at  S67  Hertz 


TABLE  7.  H-frame  Eigenvalues  for  576  Hertz 


***************** 

***  CRITERIA  TO 
***************** 


***************** 

JUDGE  RANK  *** 
***************** 


NO.* 

1 

E.VAL* 

. 785431E+02#* 

NO.* 

2 

E.VAL* 

. 393105E+02@*** 

NO.* 

3 

E.VAL* 

. 595365E+O10*************************** 

NO.* 

4 

E.VAL* 

. 119642E+00§* 

NO.- 

5 

E.VAL* 

. 764095E-01§** 

NO.  = 

6 

E.VAL* 

. 191642E-O10* 

NO.* 

7 

E.VAL* 

. 176129E-01@* 

NO.* 

8 

E.VAL* 

.117 179E-O10* 

NO.* 

9 

E.VAL* 

. 426094E-02@* 

NO.* 

10 

E.VAL* 

. 338080E-02@* 

NO.* 

11 

E.VAL* 

. 121294E-02®* 

NO.* 

12 

E.VAL* 

. 726578E-O30* 

NO.* 

13 

E.VAL* 

. 501546E-03@* 

NO.* 

14 

E.VAL* 

. 258313E-03@* 

NO.* 

15 

E.VAL* 

. 11O423E-O30* 

NO.* 

16 

E.VAL* 

. 797348E-O40* 

NO.* 

17 

E.VAL* 

. 600097 E-04@* 

NO.* 

18 

E.VAL* 

. 26O145E-O40* 

NO.* 

19 

E.VAL* 

. 211416E-04§* 

NO.* 

20 

E.VAL* 

. 166732E-04®* 
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Figure  19.  H-fr;  Tie  Third  Mode  Shape  at  576  Hertz 


TABLE  8.  H-frame  Eigenvalues  for  581  Hertz 


********************************** 

***  CRITERIA  TO  JUDGE  RANK  *** 
********************************** 


N0.= 

1 

E.VAL= 

. 119269E+03@**** 

NO.  = 

2 

E. VAL= 

. 105354E+02@* 

NO .  = 

3 

E. VAL= 

. 76534 1E+01@** ************************* 

NO .  = 

4 

E. VAL= 

. 948675E-O10* 

NO .  = 

5 

E.VAL= 

. 73O457E-O10* 

NO.  = 

6 

E. VAL= 

.  443483E-O10* 

NO .  = 

7 

E .  VAL= 

.  874408E-02@* 

NO.= 

8 

E.VAL= 

. 512513E-O20* 

NO .  = 

9 

E.VAL= 

.  387253E-O20* 

NO.  = 

10 

E . VAL= 

. 1336O9E-O20* 

NO .  = 

11 

E. VAL= 

. 7  36777E-03@* 

NO.  = 

12 

E.VAL= 

•376O61E-O30* 

NO .  = 

13 

E.VAL= 

. 21363 lE-03@* 

NO.  = 

14 

E . VAL= 

. 144119E-03@* 

NO .  = 

15 

E . VAL= 

. 122016E-03@* 

NO.  = 

16 

E. VAL= 

. 879375E-04§* 

NO.  = 

17 

E.VAL= 

. 343842E-O40* 

NO.  = 

18 

E .  VAL= 

. 232610E-04@* 

NO.  = 

19 

E .  VAL= 

. 12O455E-O40* 

NO .  = 

20 

E . VAL= 

. 903840E-05@* 
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TABLE  9.  H-frame  Eigenvalues  for  600  Hertz 


********************************** 
***  CRITERIA  TO  JUDGE  RANK  *** 
********************************** 


NO.* 

1 

E.VAL* 

. 944627E+O10***** 

e 

NO.* 

2 

E.VAL* 

. 218838E+01@* 

§ 

NO.* 

3 

E.VAL* 

. 128650E+01@***** ***************** *****§ 

NO.« 

4 

E.VAL* 

. 509607E-01@*** 

§ 

NO.* 

5 

E.VAL* 

. 171052E-016*** 

§ 

NO.* 

6 

E.VAL* 

. 628757E-O20** 

§ 

NO.* 

7 

E.VAL* 

. 297248E-O20*** 

0 

NO.* 

8 

E.VAL* 

. 984154E-O30** 

0 

NO.* 

9 

E.VAL* 

. 491233E-O30* 

@ 

NO.* 

10 

E.VAL* 

. 326437E-03@* 

§ 

NO.* 

11 

E.VAL* 

. 292134E-O30* 

@ 

NO.- 

12 

E.VAL* 

. 168954E-O30* 

0 

NO.» 

13 

E.VAL* 

. 9O895OE-O40** 

0 

NO.* 

14 

E.VAL* 

. 456719E-O40* 

0 

NO.* 

15 

E.VAL* 

. 313949E-O40* 

0 

NO.* 

16 

E.VAL* 

. 169647E-O40* 

0 

NO.- 

17 

E.VAL* 

. 1169O9E-O40* 

0 

NO.* 

18 

E.VAL* 

. 871389E-O50*** 

0 

NO.- 

19 

E.VAL* 

.344O43E-O50* 

0 

NO.* 

20 

E.VAL* 

. 244828E-O50* 

SELECTED 

RANK  OF 
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Figure  25.  H-frame  Third  Mode  Shape  at  600  Hertz 


TABLE  10.  H-frame  Modal  Assurance  Criterion  for  Transformed  Modes 


REFERENCE 

MODE 

1 

1 

1 

2 

2 

3 


ANALYSIS 

MODE 

1 

2 

3 

2 

3 

3 


M.A.C. 

1.0000000 

.0000000 

.0000000 

1.0000000 

.0000000 

1.0000000 


M.S.F. 

(REAL) 

1.0000000 

.0000000 

.0000001 

1.0000000 

-.0000002 

1.0000000 


M.S.F. 

(IMAG) 

0.0000000 

.0000000 

.0000000 

0.0000000 

0.0000000 

0.0000000 


M.S.F  Modal  Scaling  Factor  (complex) 
M.A.C.  Modal  Assurance  Criterion 


Next,  the  Multi-Mac  calculation  was  tried  on  a  circular  plate  that  is  expected  to  have  repeated  roots. 
Again.  6  reference  accelerometers  were  attached  to  the  structure  and  5  channels  around  the  peak  in 
the  quadrature  part  were  used  in  the  calculation.  Figure  26  is  a  typical  frequency  response  estimate 
and  Figure  27  is  the  summation  of  the  power  spectrum  of  the  quadrature  response.  Both  functions 
indicate  that  there  are  two  mode  in  the  frequency  range.  Tables  11  and  12  are  the  eigenvalues  from 
the  principal  component  analysis.  Both  analysis  show  two  significant  eigenvalues  indicating  two 
independent  modal  vectors.  Figures  28  and  29  are  the  transformed  mode  shapes  using  the  eigenvec¬ 
tors  from  the  principal  component  analysis.  Both  show  that  there  exists  repeated  roots  at  both  fre¬ 
quencies.  Tables  13  and  14  show  that  these  transformed  mode  shapes  are  orthogonal  with  each  other. 


TABLE  11.  Circular  Plate  Eigenvalues  for  1225  Hertz 

********************************** 


***  CRITERIA  TO  JUDGE  RANK  *** 

********************************** 

NO.-  1  E.VAL—  . 1753O7E+O50* 

0 

NO.- 

2 

E.VAL- 

. 600946E+04@***************************@ 

NO.= 

3 

E.VAL- 

. 8166O3E+O10* 

0 

NO.- 

4 

E.VAL- 

. 241708E+01®* 

0 

NO.— 

5 

E.VAL- 

.129267E+O10* 

£ 

NO.— 

6 

E.VAL- 

. 331926E+OO0* 

0 

NO.- 

7 

E.VAL- 

. 2O7313E+OO0* 

0 

NO.- 

8 

E.VAL— 

. 163597E+OO0* 

0 

NO.- 

9 

E.VAL- 

. 133O6OE+OO0* 

0 

NO.« 

10 

E.VAL- 

. 895497E-O10* 

0 

NO.— 

11 

E.VAL- 

.351897E-O10* 

0 

NO.- 

12 

E.VAL- 

.281372E-O10* 

0 

NO.- 

13 

E.VAL- 

. 244295E-O10* 

0 

NO.- 

14 

E.VAL- 

. 191932E-O10* 

0 

NO.  — 

15 

E.VAL- 

. 1OO598E-O10* 

0 

NO.- 

16 

E.VAL- 

. 729477E-O20* 

0 

NO.— 

17 

E.VAL- 

. 548224E-O20* 

0 

NO.- 

18 

E.VAL- 

. 3211O6E-O20* 

0 

NO .  — 

19 

E.VAL- 

. 14O953E-O20* 

0 

NO.* 

20 

E. VAL— 

. 1O9937E-O20* 

0 
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RANK  OF 
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Figure  28.  Circular  Plate  Mode  Shapes  at  1225  Hertz 


TABLE  12.  Circular  Plate  Eigenvalues  for  1330  Hertz 


********************************** 


***  CRITERIA  TO  JUDGE  RANK  *** 

********************************** 

NO.*  1  E.VAL*  . 220557E+05®* 

e 

NO.* 

2 

E.VAL- 

. 12 3 34 3E+O50********* **************** **a 

NO.« 

3 

E.VAL* 

.  45368OE+O10* 

@ 

NO.* 

4 

E.VAL* 

.  370407E+01§* 

% 

NO.- 

5 

E.VAL* 

. 882260E+00@* 

§ 

NO.- 

6 

E.VAL* 

. 255783E+OO0* 

§ 

NO.- 

7 

E.VAL- 

.127655E+00@* 

e 

NO.- 

8 

E.VAL- 

. 650219E-01@* 

e 

NO.« 

9 

E.VAL- 

. 469805E-01@* 

§ 

NO.- 

10 

E.VAL* 

. 386407E-01@* 

§ 

NO.- 

11 

E.VAL* 

. 225943E-01@* 

§ 

NO.- 

12 

E.VAL* 

. 209756E-01@* 

§ 

NO.- 

13 

E.VAL* 

.936288E-02@* 

e 

NO.- 

14 

E.VAL* 

. 4 8 89 2 IE-020* 

§ 

NO.- 

15 

E.VAL* 

. 391385E-02@* 

§ 

NO.« 

16 

E.VAL* 

. 283658E-02@* 

@ 

NO.- 

17 

E.VAL* 

. 189168E-02@* 

§ 

NO.- 

18 

E.VAL* 

. 150162E-02@* 

§ 

NO.- 

19 

E.VAL* 

. 123286E-02§* 

§ 

NO.- 

20 

E.VAL* 

. 593683E-03®* 

§ 

SELECTED 

RANK  OF 

DD  MATRIX*  2 

Figure  29.  Circular  Plate  Mode  Shapes  at  1330  Hertz 


TABLE  13.  Circular  Plate  Modal  Assurance  Criterion  for  Transformed  Modes  1225  Hertz  Mode 


REFERENCE 

ANALYSIS 

M.A.C. 

M.S.F. 

M.S.F. 

MODE 

MODE 

(REAL) 

(IMAG) 

1 

1 

1.0000000 

1.0000000 

1.0000000 

1 

2 

.0000000 

-.0000002 

0.0000000 

2 

2 

1.0000000 

1.0000000 

0.0000000 

M.S.F  Modal  Scaling  Factor  (complex) 
M.A.C.  Modal  Assurance  Criterion 


k 

k 

k 

k 

S 
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TABLE  14.  Circular  Plate  Modal  Assurance  Criterion  for  Transformed  Modes  1330  Hertz  Mode 


REFERENCE 

ANALYSIS 

M.A.C. 

M.S.F. 

M.S.F. 

MODE 

MODE 

(REAL) 

(IMAG) 

1 

1 

1.0000000 

1.0000000 

1.0000000 

2 

2 

1.0000000 

1.0000000 

0.0000000 

M.S.F  Modal  Scaling  Factor  (complex) 
M.A.C.  Modal  Assurance  Criterion 


5.  CONCLUSIONS 


One  of  the  conclusions  reached  in  a  previous  Air  Force  Contract  (F33615-77-C3059)  was  that  the 
area  of  modal  parameter  estimation  was,  to  in  the  future,  advance  rapidly  due  to  technology  transfer 
from  other  fields  involved  in  parameter  estimation.  This  certainly  has  occurred  as  indicated  by  the 
drastic  increase  in  the  number  of  parameter  estimation  algorithms  which  have  been  described  in  the 
literature  in  the  last  five  years.  This  effort  has  been  international  in  scope,  with  many  of  the  newer 
techniques  being  variations  of  each  other. 

These  methods  range  from  single  reference  single  degree-of-freedom  (SDOF)  methods  to  sophisti¬ 
cated  multi-reference,  multi-response,  multiple  degree-of-freedom  (MDOF)  methods.  The  algorithm 
of  choice  depends  upon  a  number  of  conditions: 


•  Modal  Application 

•  Trouble  Shooting--For  many  of  the  problem  solving,  or  trouble  shooting  applications,  the 
simplier  SDOF,  or  single  reference  MDOF  methods  are  used,  since  simple  test  procedures 
and  a  quick  look  are  desirable. 

•  Model  Verification-There  has  been  increased  emphasis  on  finite  element  verification.  These 
applications  require  a  higher  level  test  and  parameter  identification  procedures. 

•  Model  Generation,  or  Correction-There  is  also  increased  emphasis  on;  the  generation  of 
modal  models  based  upon  experimental  data,  and/or  the  correction  of  existing  models.  These 
applications  require  the  highest  level  of  test  and  parameter  identification  procedures. 

•  Equipment  Considerations 

The  parameter  identification  methods  reviewed  in  this  report  depend  heavily  upon  the  testing 
methods  (single  input,  or  multiple  input)  and  testing  equipment.  These  new  algorithms  place  a 
severe  requirement  upon  the  testing  methods  to  obtain  consistent  data  bases,  particularly  for  the 
more  advanced  multi-input  multi-output  methods. 

•  Wideband  vs  Narrowband 

Wideband  verses  narrowband  refers  to  the  frequency  bandwidth  of  the  frequency  response  meas¬ 
urements.  In  general,  for  very  broad  frequency  range  measurements,  time  domain  algorithms 
work  well,  while  frequency  domain  algorithms  seem  to  perform  well  for  the  narrow,  or  zoom 
bands.  Recently,  there  has  been  increased  emphasis  in  sine  testing.  Sine  testing,  not  in  the  classi¬ 
cal  sense,  but  in  terms  of  multi-input  multi-output  test  and  parameter  estimation  methods.  This 
emphasis  will  provide  the  impetus  to  refine  the  frequency  domain  algorithms  to  efficiently  use  the 
increased  spatial  information  that  multi-input  multi-  output  sine  testing  yields. 

•  Modal  Density 

The  choice  of  the  parameter  estimation  method  depends  heavily  upon  the  modal  density.  For 
cases  with  low  modal  density,  single  input  SDOF  or  MDOF  methods  work  well.  For  the  high 
modal  density  cases  the  multi-input  methods,  especially  ones  which  use  spatial  information,  are 
the  methods  of  choice.  It  should  be  noted  that  the  advanced  methods  require  consistent  data  and 
place  additional  constraints  on  the  testing  methods. 


A  summary  of  the  characteristics  of  the  modal  parameter  identification  methods  is  shown  in  Table 
IS.  All  of  the  methods  which  were  discussed  in  detail  in  this  report  are  briefly  summarized  in  this 
table. 


It  should  be  again  noted  that  all  of  the  methods  covered  in  this  report  can  be  described  in  terms  of  a 
characteristic  space,  where  a  particular  parameter  identification  algorithm  uses  as  input,  measured 
values  in  this  characteristic  space,  to  deconvolve  the  systems  characteristics.  The  more  advanced 
methods  use  information  from  all  three  axes  of  the  characteristic  space  simultaneously.  From  the 
measurement  standpoint,  it  is  increasingly  more  important  that  the  measured  data  be  consistent. 


TABLE  15.  Summary  of  Modal  Parameter  Estimation  Methods 


Modal  Parameter  Estimation  Characteristics 


Time, 
Frequency, 
or  Spatial 
Domain 

Single  or 
Multiple 
Degrtes-of- 
Freedom 

Olobal  Modal 
Frequencies 
and  Damping 
Factors 

Repeated  Modal 
Frequencies 
and  Damping 
Facton 

Olobal 

Modal 

Vectors 

Global 

Modal 

Participation 

Factors 

Residuals 

Quadrature  AmpHtude 

Frequency 

SDOF 

No 

No 

No 

No 

No 

Kermedy-Pancu  Ode  fit 

Frequency 

SDOF 

No 

No 

No 

No 

Yes 

SDOF  Polynomial 

Frequency 

SDOF 

Yea/No 

No 

No 

No 

No 

Non-Unear  Frequency  Domain 

Frequency 

MDOF 

No 

No 

No 

No 

Yes 

Complex  Esponenciai 

Tune 

MDOF 

No 

No 

No 

No 

No 

Least  Squares  Complex  Exponential  (LSCE) 

Thne 

MDOF 

Yes 

No 

No 

No 

No 

Ibrahim  Time  Donum  (ITD) 

Time 

MDOF 

Yes 

No 

Yes 

No 

No 

Multi-reference  Ibrahim  Time  Domain  (MITD) 

Time 

MDOF 

Yet 

Yea 

Yes/No 

NofYes 

No 

Ejprnsyatem  Realization  Algorithm  (ERA) 

Time 

MDOF 

Yes 

Yes 

Yes 

Yes 

No 

— - -  -a  a -  ■  '■< 

unno|ora  nxynmn 

Frequency 

MDOF 

Yes 

No 

No 

No 

Ye* 

Multi-reference  Orthogonal  Polynomial 

Frequency 

MDOF 

Yes 

Yes 

Yes 

Yes 

Yes 

Polyreference  Tune  Domain 

Time 

MDOF 

Yes 

Yes 

No 

Yes 

No 

Potyreference  Frequency  Domain 

Frequency 

MDOF 

Yes 

Yes 

Yes 

Yes 

Yes 

Time  Domain  Direct  Parameter  Identification 

Time 

MDOF 

Yes 

Ye* 

Yes 

Yes 

No 

Frequency  Domain  Direct  Parameter  Identification 

Frequency 

MDOF 

Yes 

Yes 

Yes 

Yes 

Yes 

MultFMAC 

Spatial 

SDOF 

No 

Yes 

Yes 

No 

No 

Multi-MAC  /  CMIF  /  Enhanced  FRF 

Spatial 

MDOF 

Yes 

Yes 

Yes 

Yes 

No 

R 
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NOMENCLATURE 
Matrix  Notation 


{-} 

braces  enclose  column  vector  expressions 

{..}T 

row  vector  expressions 

[••] 

brackets  enclose  matrix  expressions 

[■■]* 

complex  conjugate  transpose,  or  Hermitian  transpose,  of  a  matrix 

If 

transpose  of  a  matrix 

[~rl 

inverse  of  a  matrix 

[••]* 

generalized  inverse  (pseudoinverse) 

size  of  a  matrix:  q  rows,  p  columns 

r-J 

diagonal  matrix 

Operator  Notation 


A‘ 

F 

F1 

H 

H1 

In 

L 

L'x 

Re ♦ jlm 

x 

x 

7 

9 

E  AtBi 

I  si 

a_ 

at 

det[..] 


complex  conjugate 
Fourier  transform 
inverse  Fourier  transform 
Hilbert  transform 
inverse  Hilbert  transform 
natural  logrithm 
Laplace  transform 
inverse  Laplace  transform 

complex  number:  real  part  "Re",  imaginary  part  "Im" 
first  derivative  with  respect  to  time  of  dependent  variable  x 
second  derivative  with  respect  to  time  of  dependent  variable  x 
mean  value  of  y 
estimated  value  of  y 

summation  of  A,  B,  from  i  -  1  to  n 

partial  derivative  with  respect  to  independent  variable  ”t” 
determinant  of  a  matrix 
Euclidian  norm 


Roman  Alphabet 


residue  for  response  location  p,  reference  location  q,  of  mode  r 

C 

damping 

COH 

ordinary  coherence  function! 

COH* 

ordinary  coherence  function  between  any  signal  i  and  any  signal  k! 

COH “ 

conditioned  partial  coherence! 

e 

base  e  (2.71828...) 

F 

input  force 

GFF 

CFF„ 

GFF* 

[GFFX] 

GXF 

GXX 

GXX„ 

h(t) 

H(s ) 
HH 


H  t(w) 


H,(u>) 

Hs(u>) 

//*(<*>) 

[/] 

j 

K 

L 

M 

Mr 

MCOH 

H 

N< 

N0 


*/ 

Rr 

s 

t 

tk 


T 


x 


X 


z 


spectrum  of  q*  referencet 

auto  power  spectrum  of  referencet 

auto  power  spectrum  of  reference  qt 

cross  power  spectrum  of  reference  i  and  reference  kt 

reference  power  spectrum  matrix  augmented  with  the  response/reference  cross 

power  spectrum  vector  for  use  in  Gauss  elimination 

cross  power  spectrum  of  response/referencet 

auto  power  spectrum  of  response! 

auto  power  spectrum  of  response  pt 

impulse  response  functiont 

impulse  response  function  for  response  location  p,  reference  location  q  t 
transfer  function! 

frequency  response  function,  when  no  ambiguity  exist,  H  is  used  instead  of  H  (<*>)t 
frequency  response  function  for  response  location  p,  reference  location  q,  when  no 
ambiguity  exist,  H„  is  used  instead  of  (w)f 

frequency  response  function  estimate  with  noise  assumed  on  the  response,  when  no 
ambiguity  exist.  is  used  instead  of  //i(w)! 

frequency  response  function  estimate  with  noise  assumed  on  the  reference,  when  no 
ambiguity  exist,  //2  is  used  instead  of  //2(w) ! 

scaled  frequency  response  function  estimate,  when  no  ambiguity  exist,  Hs  is  used 
instead  of  //s(w) ! 

frequency  response  function  estimate  with  noise  assumed  on  both  reference  and 
response,  when  no  ambiguity  exist,  Hv  is  used  instead  of  //»(<*>)♦ 
identity  matrix 

’/T 

stiffness 

modal  participation  factor 
mass 

modal  mass  for  mode  r 
multiple  coherence  function! 
number  of  modes 
number  of  references  (inputs) 
number  of  responses  (outputs) 
output,  or  response  point  (subscript) 
input,  or  reference  point  (subscript) 
mode  number  (subscript) 
residual  inertia 
residual  flexibility 
Laplace  domain  variable 
independent  variable  of  time  (sec) 
discrete  value  of  time  (sec) 
tk  *kbt 
sample  period 

displacement  in  physical  coordinates 
response 

spectrum  ofp*  response! 

Z  domain  variable 


Greek  Alphabet 


6(0 

A/ 


Dirac  impulse  function 

discrete  interval  of  frequency  (Hertz  or  cycles/sec) 


Ar 

e 

n 

K 

TAJ 


H 


nr 

K 

{4>)r 

[*] 

{0} 

Wr 

[*] 


S 

Sr 


discrete  interval  of  sample  time  (sec) 

small  number 

noise  on  the  output 

r*  complex  eigenvalue,  or  system  pole 

Ar  -  <rr  +  M 

diagonal  matrix  of  poles  in  Laplace  domain 

noise  on  the  input 

variable  of  frequency  (rad/sec) 

imaginary  part  of  the  system  pole,  or  damped  natural  frequency,  for  mode 
(rad/sec) _ 

<4  =  n,V  WrT 

undamped  natural  frequency  (rad/sec) 


scaled  p*  response  of  normal  modal  vector  for  mode  r 
scaled  normal  modal  vector  for  mode  r 
scaled  normal  modal  vector  matrix 
scaled  eigenvector 

scaled  p *  response  of  a  complex  modal  vector  for  mode  r 
scaled  complex  modal  vector  for  mode  r 
scaled  complex  modal  vector  matrix 
variable  of  damping  (rad/sec) 

real  part  of  the  system  pole,  or  damping  factor,  for  mode  r 

damping  ratio 

damping  ratio  for  mode  r 


t 


vector  implied  by  definition  of  function 


dE 
3a  i 


(a3) 


<-> 

Equations  (a3)  and  (a4)  can  be  solved  for  the  unknown  parameters  ax  and  a0. 


'  N  ' 

E* 

N 

NY*}- 

r=l 

’  H 

r=i 

f 

(JH 

(14(1'] 

E 

H 

NY,*} 

r=i 

j* 

(a5) 


(a6) 


The  at  and  a0  values  computed  from  Eqs.  (a5)  and  (a6)  represent  the  characteristics  of  a  straight 
line  which  would  "best"  describe  the  x  and  y  sets  of  values. 

Similarly,  one  could  compute  al  and  a0  parameters  with  the  criterion  that  the  sum  of  the  square  of 
errors  in  the  x  ( e,  )  would  be  minimum.  Another  approach  could  be  the  minimization  of  the  sum  of 
the  square  of  errors  in  the  x  and  y(e,). 


Figure  A-2.  Errors  in  Least  Squares  Estimation 


The  least  squares  problem  can  be  formulated  in  matrix  notation  as: 
{Y}*[X]{A} 


where, 

>1 

*i  I 

yi 

*7  1 

{K}  = 

. 

•  1*1- 

y* 

*N  i 

A-2 


(a7) 


The  equations  presented  by  matrix  Eq.  (a7)  are  in  general  a  set  of  inconsistent  and  overdetermined 
equations.  Inconsistent,  since  it  is  not  usually  possible  to  find  {A}  that  would  satisfy  all  the  individual 
equations  of  Eq.  (a7),  and  overdetermined,  since  the  number  of  equations  is  larger  than  the  number 
of  unknowns.  The  least  square  solution  of  Eq.  (a7)  is: 

\Xf{Y}-[XriX]{A}  (a8) 

and  solving  for  unknown  vector  {A} 

{A)^[[Xr[X\'l\Xf{Y}  (a9) 

provided  that  ( \Xf[X]  )4  exists.  In  the  case  that  the  inverse  of  ( \Xf[X] )  does  not  exist,  one  could 
use  numerical  techniques  to  solve  for  vector  {A}. 

In  a  more  general  case  where  matrices  A,  X,  and  Y  are  complex  valued,  the  transpose  notation,  T, 
must  be  replaced  with  hermitian  notation,  H,  in  Eqs.  (a8)  and  (a9).  Where  the  hermitian  operator, 
H,  is  the  complex  conjugate  transpose.  Hence,  the  unknown  vector  {A}  is  given  by: 

(aio) 

Individual  equations  in  matrix  Eq.  (a7)  could  be  multiplied  by  a  weighting  factor  to  give  that  equation 
more  or  less  weight  in  the  computation.  The  weighting  factors  could  be  presented  in  form  of  a  diago¬ 
nal  (NxN)  matrix,  W.  The  diagonal  element  in  row  i  represents  the  weighting  factor  corresponding 
to  equation  i,  and  the  off  diagonal  terms  are  all  zero.  Matrix  W  is  pre-multiplied  to  both  sides  of  Eq. 
(a7): 

W\{Y)  =  [W]VO{A)  (all) 

where, 

0  .  0 
pn-  0  0 

6  d  .  wj„ 

Solving  for  vector  {A}  in  Eq.  (all): 

{A} = [x]HmBm{Y}  (ai2> 

In  general  the  relationship  of  Eq.  (al)  could  be  in  the  form  of: 

y^aN^  *aN.lxlf  l  ♦  •  •  •  +  «! x  +  a„  (al3) 

A  set  of  equations  similar  to  the  formulation  above  could  be  written  and  solved  to  obtain  the  least 
squares  estimation  of  unknown  parameters  aN  ,  as.i , ....  at ,  a0. 

The  least  squares  method  stated  above  could  easily  be  extended  to  problems  involving  more  than  one 
independent  variable.  For  example,  z  could  be  expressed  in  terms  of*  and  y : 

*  =  +  (al4) 


The  corresponding  normal  equations  for  the  least  squares  problem  of  Eq.  (al4)  are: 


LVLVWWLVUWW’J'T'A'v.'Vl.V'-VL'S.-uV.vi.V^'V'.V'.N^^ 


jj^  m  £  2|v(«a  Jr  *a1xr*  «o)]  [*yj  =  0  (al5) 

(al6) 

Jj^-£2[V(<,J-,..,.,.«<}][-i|-0  (.17) 

Equations  (al5)  -  (al7)  can  be  solved  for  the  unknown  parameters  aa  ,  ax  ,  and  a0.  The  computed 
values  «3  ,  a i ,  and  a0  represent  the  characteristics  of  a  plane  which  would  "best"  describe  the x ,  y , 
and  z  sets  of  values. 

The  above  theory  and  formulation  could  be  expanded  to  least  squares  estimation  of  a  surface  and 
eventually  to  higher  order  dimensions. 

Ai  Correlation  Coefficient 


The  "goodness"  of  the  least  squares  estimation  process  is  measured  by  the  coefficient  of  correlation 
parameter  whiclys  defined  in  terms  of  total  variation  and  explained  variation.  The  total  variation  of 

y  is  defined  as  J^jv-T)3,  which  is  the  sum  of  the  squares  of  the  deviations  of  yr  from  the  mean 

rmi 

__  ft 

value,  y.  Total  variation  consists  of  two  parts:  (1)  the  explained  variation,  £(y,-y)3;  and  (2)  the 
k  _ 

unexplained  variation,  J)(yr-yr)a.  The  terms  explained  variation  and  unexplained  variation  are  used 

rai 

to  denote  the  fact  that  the  deviations  yr-y"have  a  definite  pattern,  while,  the  deviations  yr-  yr  are  ran¬ 
dom  and  unpredictable. 


}  unexplained  variation  =  yr  -  yr 
}  explained  variation  =yr  -Y 

total  variation  »  unexp.  var.  ♦  exp.  var.  *yr  -y 


x 

Figure  A.-3.  Variations  in  Data 


Therefore,  the  coefficient  of  correlation  */a,  is  defined  as: 


(  if 

E  (9rT)2 


E  (yrJ? 


(a!8) 


The  magnitude  of  73  varies  between  0  and  1.  A  value  of  0  indicates  no  correlation  between  depen¬ 
dent  and  independent  variable(s),  while  a  value  of  1  indicates  perfect  correlation. 


It  should  be  pointed  out  that  the  coefficient  of  correlation  computed  for  a  set  of  data  and  a  assumed 
model,  only  indicates  the  relationship  of  data  based  on  the  assumed  model.  That  is,  the  coefficient  of 
correlation  measures  the  degree  to  which  the  assumed  model  describes  the  relationship  for  a  set  of 
data 


A3  Examples 


A.3.1  Example  1 

For  the  data  given  in  Table  A-l  and  the  assumed  model  equation: 


y  =  ♦ « o 


(a!9) 


X 

65 

63 

67 

64 

68 

62 

70 

66 

68 

67 

69 

71 

y_ 

68 

66 

68 

65 

69 

66 

68 

65 

71 

67 

68 

70 

TABLE  A-l.  x  and  y  Values  for  Least  Squares  Fit 


a)  Find  the  least  squares  solution  of  ax  and  a0. 

b)  Find  the  coefficient  of  correlation. 

Solution: 

a)  The  following  normal  equations  must  be  solved  for  a1  and  a0 

N  N 

«oW  +  «iE*r  =  E*  (a20) 

N  rai  N  ,=l  N 

a0E*r  +  «iE*?  =  E^r  (a21) 

r=l  rs 1  ml 

substituting  the  appropriate  terms  in  Eqs.  (a20)  and  (a21): 

12 a0  ♦  800®!  =  811 
800ao  +  53418a!  =  54107 

from  which  we  find  ax  =  0.476  and  a0  =  35.82  or 


y  =  0.476x  +  35.82 


(a22) 


yvwwn.  wv.ttuwv  vwww. vv 
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Figure  A-4.  Least  Squares  Fit  of  Data 


s 

b)  Explained  variation  »  J](yr-y  )2  =  19.22 

f=4 

s 

Total  variation  =  £(  yrY)2  -  38.92 

r=l 

10  22 

Coefficient  of  correlation  =  7a  ■  =  0.7027 

A32  Example  2 

Given  the  model  for  the  sampled  impulse  response  function  between  two  points  on  a  structure  as, 

(823) 

rsl 

where, 

*  -  0,  1,  2,  •  •  •  ,  2Af 
Af  =  value  of  time  subinterval 
U  =*Af 

and  the  known  pole  information  A,,  formulate  the  least  squares  solution  of  estimating  the  residues, 

Ar 

Solution: 

For  simplicity  and  conciseness  let, 

Zr  =eX'*‘ 

and  therefore  Eq.  (a23)  can  be  rewritten  as. 


5 


1 


a 


w 

V 

v 

s 

\ 


a 


A-6 


cvcysv;^: 


(*24) 


expanding  Eq.  (a24)  for  time  values  of  t0  to  iw.,  will  result  in  the  following  2N  equations 
h(t0)  mAx  ♦*"  +4 3^ 

ft^i)  x  Ax  Zx  ♦  A3Zj  *  ■  ■  ■  +  A 2N  Z-w 

h(t2)  -Ax  z?  *A2  z\  *Awzin 

(*25) 


ft  (tw.  i  )XAX  z?'  *  A  a  zf1'1 1  ♦  •  •  *A  w  z&J 
presenting  the  2N  equations  of  Eq.  (a25)  in  matrix  form  gives. 


1 

1 

f/lil 

ft(t  o) 

*3 

.  Z-iN 

1 

ft(t  l) 

»  *  ! 

h(t2) 

A-ui 

ft  (tin. i) 

(*26) 


or. 

Solving  Eq.  (a26)  with  more  than  2N  rows  for  vectOT  {A}  will  result  in  the  least  squares  solution  of 
the  residues. 
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APPENDIX  B:  SINGULAR  VALUE  DECOMPOSITION 
B.1  Singular  Value  Decomposition 

The  singular  value  decomposition  will  decompose  a  matrix  into  the  simplest  possible  form,  that  being 
diagonal.  This  decomposition  will  always  be  possible  regardless  of  the  rank,  or  dimension,  of  the 
matrix!1]. 

Consider  the  right  and  left  eigenvectors  of  the  m  x  n  matrix  [A  ],  which  is  of  rank  *. 


[A  ]{v} =  {u}o 

(bl) 

[A  ]*{u}*{v}<r 

(b2) 

where: 

•  {u}  -  m  x  1  left  singular  vector 

•  {  v}  =  n  x  1  right  singular  vector 

•  a  =  scalar  singular  value 

By  substituting  Eq.  (bl)  into  Eq.  (b2)  for  {«},  the  right  singular  vectors  can  be  determined  from: 

[A]*  [A]{v}*{v}<xa  (b3) 

\[A)*[A]-oni]\=0 

where: 

.1  =  1-.*  of>0 
.  i  =*+l  -*n  a f  -  0 

•  {vfo  {v}3  {v}„]  right  singular  unitary  matrix. 

By  substituting  Eq.  (b2)  into  Eq.  (bl)  for  { v},  the  left  singular  vectors  can  be  determined  from: 

[A]  [A  ]*{«}*{«}  *a  0>4) 

|  [A  ]  [A  ]*  -<r3[7]  |  =0 

where: 

.i  =  l-*  of>  0 

•  i  =*+l—  m  of  =0 

•  =  {«}a  {«}■,  ]  left  singular  unitary  matrix. 

If  the  eigenvector  matrices  [  t/  ]  and  [  V ]  are  unitary  ([£/]*  m  [t/f1  and  [  K  ]ff  *  [  K  J"1 ),  that  is,  both 
columns  and  rows  form  an  orthonormal  set,  the  linear  transformation  will  preserve  both  angles  and 
lengths.  Interpreted  geometrically,  linear  transformations  defined  by  unitary  matrices  behave  like 
simple  rotations  in  space.  The  matrix  [A  ],  can  now  be  decomposed  into  diagonal  form  by  appending 
the  eigenvector  matrices  to  Eq.  (bl)  and  premultiplying  by  [  U\B,  forming  the  matrix  product: 
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i  “ 
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B-l 


(b5) 
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•••  0 

0 

<r3 
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•••  0  [0] 
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0 

^3 

0 

6 

6 

0 

•  Ok 

— 

m 

— 

“  ■  m 

-  {“}f  [{«h*i  ••  {«}***  {0}*+i  •  •  {0},]  . 

{«}£ 

Using  the  unitary  matrix  property,  [  U  ]B  [  U  ]-[/],  the  right  hand  side  can  be  further  simplified  as: 


tra  foil  ■ 


Thus,  the  origonal  matrix,  [A  ],  is  decomposed  into  the  matrix  [  S  ],  with  the  singular  values  on  the 
diagonal,  by  the  unitary  matrices  [  U  ]  and  [  V]. 

Then  the  singular  value  decomposition  of  [A  ]  is  defined  by: 

Noting  the  partitioning  of  the  matrix  [5  ]  of  Eq.  (b6),  the  unitary  transformation  matrices  [  U ]  and 
[  V]  can  be  partitioned  in  the  same  way  to  yield: 

[tfla][[^  [jjj]  [yfa  (b8) 

Thus,  the  singular  value  decomposition  of  [A  ]  can  be  further  reduced  to: 

where: 

•  [U]j  =  left  singular  submatrix  of  size  mxk 

•  f  EJ  =  diagonal  singular  value  matrix  of  size  kxk 

•  [ v  Jf  =  right  singular  submatrix  of  size  kxn. 

The  (Moore-Penrose)  generalized  inverse  (pseudoinverse)  l1'3!,  [A  J\  of  [A  ]  is: 

[oj  mB=[vw\-'[u]f . 

where: 

•  [A  }*  «  n  xm  generalized  inverse  of  [A  ] 

•  [y] j  *  right  singular  submatrix  of  size  nxk. 

•  f  EJ'1  »  inverse  of  diagonal  singular  value  matrix  of  size  kxk 

•  l  U  ]f  *  left  singular  submatrix  of  size  kxm. 
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